Continuum Limit Physics from 2+1 Flavor Domain Wall QCD
Y.Aoki, R.Arthur, T.Blum, P.A.Boyle, D.Brommel, N.H.Christ, C.Dawson, J.M.Flynn, T.Izubuchi, X-Y.Jin, C.Jung, C.Kelly, M.Li, A.Lichtl, M.Lightman, M.F.Lin, R.D.Mawhinney, C.M.Maynard, S.Ohta, B.J.Pendleton, C.T.Sachrajda, E.E.Scholz, A.Soni, J.Wennekers, J.M.Zanotti, R.Zhou
aa r X i v : . [ h e p - l a t ] J u l CU-TP-1195, Edinburgh 2010/13, KEK-TH-1380, RBRC-847, SHEP-1027
Continuum Limit Physics from 2+1 Flavor Domain Wall QCD
Y. Aoki, R. Arthur, T. Blum,
3, 1
P.A. Boyle, D. Br¨ommel,
4, 5
N.H. Christ, C. Dawson, J.M. Flynn, T. Izubuchi,
1, 8
X-Y. Jin, C. Jung, C. Kelly, M. Li, A. Lichtl, M. Lightman, M.F. Lin, R.D. Mawhinney, C.M. Maynard, S. Ohta,
11, 12, 1
B.J. Pendleton, C.T. Sachrajda, E.E. Scholz, A. Soni, J. Wennekers, J.M. Zanotti, and R. Zhou (RBC and UKQCD Collaborations) RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA SUPA, School of Physics, The University of Edinburgh, Edinburgh EH9 3JZ, UK Physics Department, University of Connecticut, Storrs, CT 06269-3046, USA School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK J¨ulich Supercomputing Centre, Institute for Advanced Simulation,Forschungszentrum J¨ulich GmbH, 52425 J¨ulich, Germany Physics Department, Columbia University, New York, NY 10027, USA Dept. of Physics, University of Virginia, 382 McCormick Rd. Charlottesville, VA 22904-4714 Brookhaven National Laboratory, Upton, NY 11973, USA Department of Physics, Yale University, Sloane Physics Laboratory, New Haven, CT 06511, USA EPCC, School of Physics, The University of Edinburgh, Edinburgh EH9 3JZ, UK Institute of Particle and Nuclear Studies, KEK, Tsukuba, Ibaraki, 305-0801, Japan Department of Particle and Nuclear Physics, Sokendai GraduateUniversity of Advanced Studies, Hayama, Kanagawa 240-0193, Japan Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany
PACS numbers: 11.15.Ha, 11.30.Rd, 12.15.Ff, 12.38.Gc 12.39.Fe
ABSTRACTWe present physical results obtained from simulations using 2+1 flavors of domain wall quarksand the Iwasaki gauge action at two values of the lattice spacing a , ( a − = 1.73 (3) GeV and a − = 2.28 (3) GeV). On the coarser lattice, with 24 × ×
16 points (where the 16 correspondsto L s , the extent of the 5 th dimension inherent in the domain wall fermion (DWF) formulationof QCD), the analysis of ref. [1] is extended to approximately twice the number of configura-tions. The ensembles on the finer 32 × ×
16 lattice are new. We explain in detail how we uselattice data obtained at several values of the lattice spacing and for a range of quark masses incombined continuum-chiral fits in order to obtain results in the continuum limit and at physicalquark masses. We implement this procedure for our data at two lattice spacings and with uni-tary pion masses in the approximate range 290–420 MeV (225–420 MeV for partially quenchedpions). We use the masses of the p and K mesons and the W baryon to determine the physicalquark masses and the values of the lattice spacing. While our data in the mass ranges above areconsistent with the predictions of next-to-leading order SU(2) chiral perturbation theory, they arealso consistent with a simple analytic ansatz leading to an inherent uncertainty in how best to per-form the chiral extrapolation that we are reluctant to reduce with model-dependent assumptionsabout higher order corrections. In some cases, particularly for f p , the pion leptonic decay con-stant, the uncertainty in the chiral extrapolation dominates the systematic error. Our main resultsinclude f p = ( ) stat ( ) syst MeV, f K / f p = . ( )( ) where f K is the kaon decay constant, m MS s ( ) = ( . ± . ) MeV and m MS ud ( ) = ( . ± . ) MeV ( m s / m ud = . ± . m s and m ud are the mass of the strange-quark and the average of the up and down quarkmasses respectively, [ S MS ( )] / = ( ) MeV, where S is the chiral condensate, the Som-mer scale r = . ( ) fm and r = . ( ) fm. I. INTRODUCTION
For several years now, the RBC and UKQCD Collaborations have been undertaking a major pro-gramme of research in particle physics using lattice QCD with Domain Wall Fermions (DWF)and the Iwasaki gauge action. In the series of papers [1–3], we studied general properties of en-sembles with an inverse lattice spacing of a − = . ( ) GeV (corresponding to b = .
13) andwith unitary pion masses m p ≥
330 MeV (partially quenched m p &
240 MeV). The number ofpoints in these ensembles are 16 × × × ×
16 [3] and 24 × ×
16 [1], where thefifth dimension is a feature of DWF and is not visible to low-energy physics which remains four-dimensional. We do not review the properties of DWF here, beyond underlining their physicalchiral and flavor properties which we exploit in much of our wider scientific programme. We haveused these ensembles to investigate a broad range of physics, including studies of the hadronicspectrum, mesonic decay constants and light-quark masses [1], the evaluation of the B K parameterof neutral-kaon mixing [1, 4], the calculation of the form-factors of K ℓ decays [5, 6], studies innucleon structure [7–9] and proton decay matrix elements [10] and very recently the first latticestudy of the masses and mixing of the h and h ′ mesons [11] as well as a determination of thematrix elements relevant for neutral B -meson mixing in the static limit [12]. A key limiting factorin the precision of these results was that the simulations were performed at a single lattice spacing.In this paper we remove this limitation, by presenting results for the spectrum, decay constants andquark masses obtained with the same lattice action using ensembles generated on a 32 × × b = .
25, for which we will seebelow that a − = . ( ) GeV. Now that we have results for the same physical quantities with thesame action at two values of the lattice spacing we are able to perform a continuum extrapolationand below we will present physical results in the continuum limit.Since the most precise results at b = .
13 were obtained on the 24 × ×
16 [1] lattices, as ashorthand throughout this paper we will refer to these lattices as the 24 ensembles and label thenew lattices at b = .
25 as the 32 ensembles.The new 32 ensembles at b = .
25 will, of course, be widely used also in our studies of otherphysical quantities. In this first paper however, we discuss their properties in some detail (seeSec. II). In this section we also discuss reweighting which allows us to eliminate one source ofsystematic uncertainty. While at present we cannot simulate with physical u and d quark masses,there is no reason, in principle, why we cannot simulate with the physical strange quark mass.The difficulty however, is that we don’t know a priori what this mass is and so in practice thesimulations are performed with a strange quark mass which is a little different from the physicalone. As explained in Section II D, the technique of reweighting allows us to correct a posteriorifor the small difference in the simulated and physical strange quark masses. In Section III, wepresent updated raw results for the pion and kaon masses and decay constants and the mass ofthe W -baryon on the 24 ensembles which have been extended beyond those discussed in ref.[1].Section IV contains the corresponding results on the 32 ensembles. In these two sections we alsopresent the raw results for the masses of the nucleon and D baryons from the two ensembles, but incontrast to the mesonic quantities a description of their chiral behaviour and extrapolation to thecontinuum limit are postponed to a future paper.The price we pay for using a formulation with good chiral and flavor properties is the presenceof the fifth dimension and the corresponding increase in computational cost. The lightest unitarypion which we have been able to afford to simulate has a mass of 290 MeV and so, in additionto the continuum extrapolation we need to perform the chiral extrapolation in the quark masses.In Sec. V we present a detailed explanation of how we combine the chiral and continuum extrap-olations in an attempt to optimize the precision of the results, exploiting the Symanzik effectivetheory approach as well as chiral perturbation theory and other ansatze for the mass dependence ofphysical quantities. Having explained the procedure, we then proceed in Section V E to discuss theresults, to determine the physical bare masses and lattice spacings as well as to make predictionsfor the pion and kaon decay constants. In particular we find that the ratio of kaon and pion decayconstants [78] f K f p = . ± . , (1)where the error is largely due to the uncertainty in the chiral behaviour of f p as explained inSec. V E 3. From the chiral behaviour of the masses and decay constants we determine the corre-sponding Low Energy Constants (LECs) of SU(2) Chiral Perturbation Theory (ChPT).Among the most important results of this paper are those for the average u and d quark mass andfor the strange quark mass which are obtained in Sec.VI: m MS ud ( ) = ( . ± . ) MeV and m MS s ( ) = ( . ± . ) MeV . (2)The masses are presented in the MS scheme at a renormalization scale of 2 GeV, after the renor-malization to symmetric momentum schemes has been performed non-perturbatively [13, 14] andthe conversion to the MS scheme has been done using very recent two-loop results [15, 16].Section VII contains a discussion of the topological charge and susceptibility of both the 24 and32 ensembles and in Sec.VIII we summarise our main results and present our conclusions. Thereare three appendices. Appendix A contains the chiral extrapolations performed separately onthe 24 and 32 ensembles. This is in contrast with the procedure described in Section V E inwhich the chiral and continuum extrapolations were performed simultaneously with common fitparameters at the two spacings. Appendix B contains a detailed analysis of a subtle issue, thenormalization of the partially conserved axial current. For domain wall fermions this is expectedto deviate from the conventionally normalized continuum current by terms of order am res , where a is the lattice spacing and m res is the residual mass [1, 17]. Current simulations are now becomingsufficiently precise that these effects need to be understood and quantified and the method proposedin appendix B, in which the O ( am res ) effects are absent, is implemented in the numerical analysesthroughout the paper. Finally Appendix C contains a discussion of the expected statistical errorswhen reweighting is performed on Monte Carlo data to obtain results with a different action fromthat used to generate the data.We end the Introduction with an explanation of our notation for quark masses [1]. When discussingunitary computations, with the valence and sea quarks degenerate, we call the bare light ( u or d )quark mass m l and the bare heavy (strange) quark mass m h . m ud and m s refer to the physical valuesof these masses (we work in the isospin limit so that the up and down quarks are degenerate). Forthe partially quenched computations we retain the notation m l and m h for the sea-quark masses,but use m x and m y for the valence quarks. A tilde over the mass indicates that the residual mass has been added, e m q = m q + m res ; it is e m which is multiplicatively renormalizable. II. SIMULATION DETAILS AND ENSEMBLE PROPERTIES
As described in Ref. [1, 3, 18], we generate ensembles using a combination of the DWF formula-tion of Shamir [19] and the Iwasaki gauge action [20]. For the fermionic action we use a value of1.8 for the “domain wall height” M and an extension of the 5 th dimension of L s =
16. In additionto the new ensembles generated on a 32 ×
64 lattice volume and a gauge coupling b = . × b = .
13 ensembles generated in our previousstudy [1]. As indicated in Tab. I we have extended the m l = . ×
64 ensemble from 4460to 8980 MD units while the m l = .
01 ensemble has been extended from 5020 to 8540 MD units.The three 32 ×
64 ensembles that are first reported here are also shown in Tab. I and those withlight quark masses of 0.004, 0.006 and 0.008 contain 6856, 7650 and 5930 MD units respectively.
A. Ensemble Generation
For the generation of both the 24 ×
64 and 32 ×
64 ensembles, we employ the “RHMC II”algorithm described in Ref. [1]. More specifically, the simulation of two light quarks and onestrange quark is carried out using a product of three separate strange quark determinants eachevaluated using the rational approximation. The 2 flavors of light quarks are preconditioned bythe strange quark determinant [21]. While the preconditioning mass does not have to be the sameas the strange-quark mass, we found that the strange-quark mass is close to being optimal in DWFsimulations in tests on smaller volumes.Using the notation D ( m l ) = D † DW F ( M , m l ) D DW F ( M , m l ) , the fermion determinant including thecontribution from the Pauli-Villars fields and evaluated on a fixed gauge configuration can bewritten asdet " D ( m s ) / D ( m l ) D ( ) / = det (cid:20) D ( m s ) D ( ) (cid:21) / · det (cid:20) D ( m l ) D ( m s ) (cid:21) (3) = det (cid:20) R (cid:18) D ( m s ) D ( ) (cid:19)(cid:21) · det (cid:20) R (cid:18) D ( m s ) D ( ) (cid:19)(cid:21) · det (cid:20) R (cid:18) D ( m s ) D ( ) (cid:19)(cid:21) · det (cid:20) D ( m l ) D ( m s ) (cid:21) . (4)In the third line we explicitly show how this ratio of determinants is implemented using the ratio-nal approximation. Here R a ( x ) denotes x a evaluated using the rational approximation and eachdeterminant is evaluated using a separate set of pseudofermion fields. An Omelyan integrator [22]with the Omelyan parameter l = .
22 was used in each part of evolution.Given the disparate contributions to the molecular dynamics force coming from the gauge actionand the different factors in Eq. (4) we follow the strategy of Ref. [23] and increase performance bysimulating these different contributions with different molecular dynamics time step granularities.In particular, the suppression of the force from the light quark determinant that results from theHasenbusch preconditioning allows us to evaluate the computationally expensive force from thelight quark using the largest time step among the different terms, decreasing the computational costsignificantly. As a result, we divide our simulation in such a way that D t light : D t heavy : D t gauge = / D t heavy effectively halfof D t light .) This ratio of time steps was used for all the ensembles studied here. However D t light was varied from ensemble to ensemble to reach an approximate acceptance of 70%. The precisenumbers that were used are listed in Tab. I.In addition, we chose to simulate with a trajectory length t = ensembles, twicethat used for the 24 ensembles. While a longer trajectory length may be expected to reduce theautocorrelation between configurations, the time for a trajectory scales very nearly linearly in thetrajectory length. In comparisons between t = t = m s a m l a ˜ m s / ˜ m l D t light t (Ref.[1]) t (MD) Acceptance h P i h ¯ yy ( m l ) i V / a = × L s = , b = . , a − = . ( ) GeV, m res a = . ( ) , t / traj = V / a = × L s = , b = . , a − = . ( ) GeV, m res a = . ( ) , t / traj = h P i ) and value for the light-quark chiral condensate ( h ¯ yy ( m l ) i ) for the ensembles studied in this paper. The fifth column shows thenumber of time units in the ensembles that were included from Ref. [1]. The residual masses given explicitlyand those appearing in the ratio e m l / e m s are taken from Table VII appearing in Section III below. A final optimization was used for the simulations run on the IBM BG/P machines at the Ar-gonne Leadership Computing Facility(ALCF). Instead of using double precision throughout, theBAGEL-generated assembly routines [24] keep the spin-projected spinors in single precision inthe conjugate gradient(CG) inverters during the molecular dynamics evolution to decrease theamount of communication needed per CG iteration. (Full precision is used in the accept-rejectstep.) While this kind of improvement is expected to make the molecular dynamics integrator un-stable for sufficiently large volumes, the effect on the acceptance turned out to be minimal for allthe ensembles presented in this paper while improving the performance of the CG by up to 20%compared to a full double precision CG with the same local volume. l = 0.0040.6150.61550.616 m l = 0.0060 2000 4000 6000 8000MD units0.6150.61550.616 m l = 0.008 0.00050.00060.00070.0008 m l = 0.0040.00070.00080.00090.001 m l = 0.0060 1000 2000 3000 4000 5000 6000 7000 80000.00090.0010.00110.0012 m l = 0.008 FIG. 1: Evolution of the average plaquette (left panel) and the chiral condensate (right panel) for the b = . , × L s =
16 ensembles. The chiral condensate is normalized such that h ¯ yy i ∼ B. Ensemble properties
In Fig. 1 we show the evolution of the plaquette and the chiral condensate for the 32 ensembles.Both quantities suggest that 500 MD units is enough for the thermalization of each of the 32 ensembles. We have thus begun measurements at 1000 MD units for m l = .
006 (except for themeasurements of the chiral condensate which started after 3304 MD units) and 520 MD units forthe other 32 ensembles. (The starting points for measurements on the three 24 ×
64 ensemblesare given in Tab. I of Ref. [1].)Figure 2 shows the integrated autocorrelation time for various quantities measured on the 32 ensembles. As can be seen the plaquette, chiral condensate and even the light pion propagatorfor a separation of 20 time units show a short autocorrelation time of 5-10 MD units. However,the measured autocorrelation times for the topological charge are much larger, on the order of80 MD units. In fact, as is discussed in Section VII, the evolutions shown in Fig. 52 suggesteven longer autocorrelation times implying that the autocorrelation times shown in Fig. 2 may beunderestimated because of insufficient statistics.In Section VII this issue of the autocorrelation time for the topological charge is discussed ingreater detail and the b = .
13 and 2.25 evolutions are compared. The 32 , b = .
25 ensem-bles (with finer lattice spacing) are shown to evolve topology more slowly. This suggests that
Plaquette Chiral condensate Pseudoscalar(t=20) Topological charge(m l =0.004) Topological charge(m l =0.006) Topological charge(m l =0.008) FIG. 2: The integrated autocorrelation time is shown for the average plaquette, chiral condensate h ¯ yy i ,pseudoscalar propagator at time separation 20 from a Gaussian source and point sink, all computed fromthe 32 , m l = .
004 ensemble and the global topological charge for all three 32 ensembles. The chiralcondensate and plaquette are measured every two MD units and the averages within sequential blocks of10 MD units have been analyzed. The topological charge is measured every 4 MD units and the averageswithin sequential blocks of 20 MD units have been analyzed. All other quantities were measured every 20MD units and no averaging has been performed. Further discussion of the topological charge is given insection VII. the change from the DBW2 gauge action used in earlier 2-flavor work [25] to the Iwasaki gaugeaction used here may have been a wise one. While the DBW2 gauge action gives smaller residualDWF chiral symmetry breaking, it does this by suppressing the tunneling which changes topolog-ical charge. Thus, the use of the DBW2 gauge action may have resulted in a topological chargeevolution for our current finest lattice spacings that would have been unacceptably slow. C. Fitting procedure
In the analysis described in this paper it is important to take into account the fact that the variousquantities computed on a single gauge configuration may be correlated. To do this we apply thejackknife technique to simple uncorrelated fits. While there is no proof, or even expectation, thatthis is an optimal procedure, the jackknife will provide a good estimate of the error except in the0unlikely event of large deviations of our result from a normal distribution. While we could attemptto perform a “text-book” correlated fit (again, using a jackknife procedure), this would not besensible: such fits assume that the data should exactly follow the functional form used in the fit.In the case of a fit to chiral perturbation theory or a simpler analytic ansatz for the quark-massdependence of physical quantities we know that this is not the case. While this complaint appliesto both correlated and uncorrelated fits, for the highly correlated lattice data with which we aredealing, small deviations (which in this procedure are assumed to be statistical, but in our case arelikely to be systematic) are penalized by many orders of magnitude more for the correlated thanuncorrelated fits. Nevertheless, we have performed correlated fits, where the correlation matrixis obtained by taking increasing numbers of the leading eigenvectors. Within our limited abilityto estimate the correlation matrix, we find no significant difference in the results and errors withthose obtained using uncorrelated fits. Therefore, in this paper (as was also the case in Ref. [1])we present our main results from the uncorrelated fits, but with a full jackknife procedure forestimating the errors. However, it must be borne in mind that for such uncorrelated fits the resulting c may not be a reliable indicator of goodness of fit. Therefore, we present a sample set of our fitsgraphically. D. Reweighting in the mass of the sea strange-quark
The sea strange quark mass value used in our ensemble generation, m (sim) h , differs from the onein nature, which we determine only after performing our final analysis. In this subsection, wedescribe the reweighting method used to correct this strange quark mass from m (sim) h to the tar-get mass m h . Various target heavy quark masses are determined in Section V through interpola-tion/extrapolation to yield meson masses which match either unphysical values present in a dif-ferent ensemble or which reproduce those from experiment. Recently, several large-scale QCDsimulations have been reported using a reweighting technique [26–28]. The various uses of thismethod include obtaining sea quark mass derivatives in Ref. [29], tuning the light and strangequark masses in Ref. [30], tuning the strange and charm quark masses in Ref. [31] and going tolarger L s for the DWF action in Ref. [32].An observable, such as the meson propagator, at the target strange sea quark mass m h is obtained bymeasuring that observable on the ensemble generated using m (sim) h , multiplied by the reweighting1factor w : h O i m h = h O w i m (sim) h h w i m (sim) h . (5)Here the reweighting factor w [ U m ] for a particular ensemble of gauge links U m is the ratio of thesquare root of the two-flavor Dirac determinant evaluated at the mass m h divided by that samerooted determinant evaluated at m (sim) h , w [ U m ] = det D ( m h ) / det D ( m (sim) h ) / . (6)This factor must be calculated for each configuration on which measurements will be performedin the ensemble generated using the sea strange mass m (sim) h .Among the many possible ways of computing the determinant ratio in Eq. (6), we have chosen touse the Hermitian matrix W ( m h , m (sim) h ) , whose determinant is w [ U m ] , W ( m h , m (sim) h ) = h D ( m (sim) h ) † i / h D ( m h ) † i − / [ D ( m h )] − / h D ( m (sim) h ) i / . (7)The square root of these matrices is implemented using the same rational polynomial approx-imation, R ( x ) , and multi-shift conjugate gradient algorithm, which are used in the ensemblegeneration. The order of the matrix products in W assures that in the limit of m h → m (sim) h , W goesto the unit matrix, so that the method described below for evaluating w has vanishing stochasticerror in this limit.To obtain w on each configuration, the determinant of W is stochastically evaluated using a com-plex random Gaussian vector x of dimension L s ×
12. Each complex element is drawn from arandom distribution centered at zero with width s x in both the real and imaginary directions: w = hh e − x † [ W − / ( s x )] x ii x ≡ R D x D x † e − x † [ W − / ( s x )] x e − x † x / ( s x ) R D x D x † e − x † x / ( s x ) . (8)We set s x = / N x Gaussian vectors per configuration. For one sample, twomulti-mass inversions, one for m h and another for m (sim) h , are performed.One needs to be careful in evaluating Eq. (8) to avoid a large and difficult to estimate statisticalerror. When the eigenvalues of W , l W , are far from 1 / ( s x ) , the large shift in the width of theGaussian in the integrand will cause poor sampling in this stochastic evaluation of w , as can beseen if Eq. 8 is rewritten with W diagonal: w = (cid:213) l W ∈ spect ( W ) Z d x l x † l e − x † l [ l W − / ( s x )] x l e − x † l x l / ( s x ) / (cid:213) l W ∈ spect ( W ) Z d x l x † l e − x † l x l / ( s x ) . (9)2The first exponential function in the integrand (9) will be a rapidly decreasing function of x † x when [ l W − / ( s x )] is large, with most of the Gaussian samples generated according to thesecond exponential function in Eq. (9) falling in a region where the first factor is very small. Inthis sense, Eq. (8) may provide a statistically noisy estimate of the ratio of the determinants inEq. (6). The fluctuations in this estimate will be rapidly reduced when [ l W − / ( s x )] → s x , when W becomes close to the unit matrix, W → W is divided into N rw factors [27] w = det W = N rw − (cid:213) i = det W i = N rw − (cid:213) i = hh e − x † i [ W i − / ( s x )] x i ii x i . (10)Each of W i needs to be close to unit matrix while keeping the determinant of the product thesame as the original determinant. Each factor det W i in the product, is evaluated using Eq. (8)with N x Gaussian vectors. We note that all Gaussian vectors, x i , must be statistically independent otherwise there will be unwanted correlation among contribution from the N rw steps. A similardecomposition of the reweighting factor is also possible by using the n th root of the operators[32].In this work, W i is chosen by uniformly dividing the interval [ m h , m (sim) h ] into smaller pieces: W i = W (cid:16) m ( i + ) h , m ( i ) h (cid:17) , (11) m ( i ) h = m (sim) h + i m h − m (sim) h N rw , ( i = , , · · · , N rw ) . (12)In that way, reweighting factors for the intermediate masses m ( i ) h are also obtained, which will beused in our analysis too.For a given difference between the target and the simulation masses, m h − m (sim) h , N rw needs to besufficiently large that W i is close to the unit matrix, suppressing the statistical noise in estimatingeach of the determinants. We have checked whether N rw is large enough in our calculation ofthe reweighting factor. Figure 3 shows the logarithm of the full reweighting factor, − ln ( w ) , as afunction of the number of divisions in strange quark mass, N rw , on the b = . , × , m l = .
005 lattices, the 2,000th trajectory in the left panel and the 4,000th trajectory in the right panel.The target and simulation quark masses are m h = .
035 and m (sim) h = . N rw ≤
10, the reweighting factor w appears inconsistent with the results obtained for larger N rm by a large amount (note that − ln ( w ) is plotted) for the left case (2,000th trajectory). We believethis is caused by the poor stochastic sampling in our method to compute w when N rw ≤
10 andthat for these cases the statistics are insufficient to estimate the error accurately.3 N rw - l n ( w ) traj=2000 N rw traj=4000 FIG. 3: Logarithm of the reweighting factor, − ln ( w ) , as a function of the number of divisions in thestrange quark mass, N rw on the b = . , × , m l = .
005 lattices, the 2,000th trajectory on the leftpanel and the 4,000th trajectory on the right panel. The target and simulation quark masses are m h = . m (sim) h = . N rw = , , , , ,
40, the number of Gaussian samples per mass steps is set to N x = , , , , ,
2, respectively. The error bars shown are the standard deviations resulting from N rw × N x samples for det W i . We interpret the inconsistency between the values for N rw =
1, 5 and 10 and those withlarger N rw in the left-hand panel as resulting from insufficient statistics leading to under-estimated errors forthese three cases where the stochastic sampling is very poor.ensemble m (sim) h m h N rw N x ×
64 0.030 0.025 10 424 ×
64 0.040 0.030 40 2TABLE II: Parameters chosen for the sea strange quark mass reweighting calculation.
We also check the relative difference between the reweighting factors for N rw =
20 and N rw = N rw =
20 is sufficient to estimate the reweightingfactor and its error for changing from m (sim) h = .
040 to m h = .
035 on this ensemble. We summa-rize the values of N rw and N x used in estimating the reweighting factors for the sea strange quarkmass in Tab. II.Is the N rw dependence, described above, all one needs to check to assure the correctness of the4 ( w ( )- w ( ) ) / w ( ) FIG. 4: The relative differences between the reweighting factors for N rw = , N x = N rw = , N x = m h = .
035 and m (sim) h = . reweighting procedure? The answer is clearly no. So far, we have only established that Eq. (10)estimates w to some degree of accuracy, on each configuration for large N rw . One needs furtherchecks to see whether or not the reweighted observable in Eq. (5) has an accurately estimatedstatistical error. A highly inaccurate estimate of the statistical errors could easily result from apoor overlap between the reweighted ensemble and the original ensemble generated by the RHMCsimulation. In addition, because the reweighted observable in Eq. (5) is given by a ratio of averagesit is a biased estimator of the observable of interest. In this circumstance, a large statistical error,even if well determined, may lead to a systematic error of order 1 / N conf enhanced by this largestatistical error.We have attempted the following checks: In Fig. 5, w is plotted as a function of trajectory. If thefluctuation among different configurations is large, Eq. (5) might be dominated by a small numberof measurements made on those configurations with large w , and the measurement efficiency forthe reweighted observable would be very poor. Using the reweighting factor, w i , obtained on the i th configuration, the reweighted observable O can be written from Eq. (5) as, h O i m s = N conf (cid:229) i = O i ˆ w i , (13)ˆ w i = w i (cid:229) N conf i = w i . (14)Because the process of reweighing selectively samples the original distribution, even with pre-5 ensemble max ( w i ) min ( w i ) N Eff N ∗ Eff N conf × , m l = .
005 10.0 0.078 90.3 20.3 20324 × , m l = .
010 5.50 0.049 97.0 32.4 17832 × , m l = .
004 4.77 0.17 228 63.9 30532 × , m l = .
006 3.45 0.23 234 90.4 31232 × , m l = .
008 5.36 0.16 183 47.0 252TABLE III: The maximum and minimum reweighting factors, the effective number of samples, N eff , ac-cording to the formula derived in this paper, (Eq. (15)), the corresponding number, N ∗ eff given by theformula of Ref. [33] (defined in Eq. (16)) and the actual number of configurations N conf in each en-semble. The target sea strange quark mass and that of the simulation are m h = . , m (sim) h = . m h = . , m (sim) h = . ×
64 (32 × cisely determined reweighting factors we should expect the effective number of samples to bereduced and the statistical errors to increase. In Appendix C this effect is analyzed in the case thatcorrelations between the data and the reweighting factors can be neglected when estimating thesestatistical errors, including the effects of autocorrelations. For the case of no autocorrelations, weobtain the following expression for the effective number of configurations after reweighting: N eff = (cid:16) (cid:229) N conf n = w n (cid:17) (cid:229) N conf n = w n . (15)The quantity N eff goes to N conf if there is no fluctuation in the w i while it goes to 1 if the largest w i completely dominates the reweighted ensemble. We summarize the statistical features of thereweighting factors for each ensemble in Tab. III. For completeness we also compare the definitionof N eff given in Eq. (15) with the more pessimistic estimate used in Ref. [33]: N ∗ eff = (cid:229) N conf i = w i max j ( w j ) . (16)As can be seen from Tab. III, our choice gives a somewhat more optimistic view of the effects ofreweighting on the effective size of our ensembles.As the numbers in Tab. III indicate, for our ensemble and reweighting settings, the ensembles arenot overwhelmed by a small number of configurations.The efficiency of the reweighting procedure is also observable dependent. It is influenced by thefluctuations of the reweighted observable within the ensemble and the strength of the correlation6 w norm at m h =0.0345 from 0.040 m l =0.005, 0.010 w norm at m h =0.0275 from 0.030 m l =0.004, 0.006, 0.008 FIG. 5: The normalized reweighting factor ˆ w i as a function of trajectory number i for the 24 × , m l = . , .
010 ensembles (left-hand plot) and the 32 × , m l = . , . , .
008 ensembles (right-handplot). The sea quark masses m l are plotted in ascending order from top to bottom. The target sea strangequark mass and that of simulation are m h = . , m (sim) h = .
040 ( m h = . , m (sim) h = . between the reweighted observable and the reweighting factor. Sanity checks of the statisticalproperties of the most important observables, m p and f p , have been performed and are summarizedin Fig. 6. The observables reweighted to m h = . m (sim) h = .
030 are calculated usingthe first half and the second half of the ensemble (circle symbols), which are compared to thatof the full statistics (square symbols). The number of the Gaussian vectors, N x , is also variedfrom N x = N x = m p , all thestatistical samples are within 1 × s , while for f p the deviations are less than ∼ × s .To probe the m h dependence of the observables, we show in Fig. 7 the correctly reweighted m p and f p as a function of m h along with the results obtained from randomly permuting the { w i } inEq. (13). The random permutation is done for each reweighted mass m h to show the differencefrom the correctly reweighted observables. While the randomly reweighted observables are almostflat in m h , the correctly reweighted observables have a positive slope in m h . Finally in Fig. 8 weplot the reweighted observables f p and f K as a function of the target reweighted mass m h for three7 Unitary m p versus hit index reweighted to 0.02500 hit index m p ( G e V ) single hit per conftwo hits per conffour hits per confunreweighted value 0 2 4 6 8 100.1320.1340.1360.1380.140.1420.1440.146 Unitary f p versus hit index reweighted to 0.02500 hit index f p ( G e V ) single hit per conftwo hits per conffour hits per confunreweighted value FIG. 6: Reweighted values for m p (left) and f p (right) for various numbers of reweighting hits, N x = N x = N x = × × , ( m l , m h ) = ( . , . ) ensemble with a light valence quark of mass 0.004. Theblack symbols are the unreweighted observables. Unitary m π versus m h from reweighting m h total (GeV) m π ( G e V ) Randomly permuted reweighing factorsCorrect reweighting factors 0.09 0.095 0.1 0.105 0.11 0.1150.1350.140.145
Unitary f π versus m h from reweighting, 0.004/0.03 m h total (GeV) f π ( G e V ) Randomly permuted reweighing factorsCorrect reweighting factors
FIG. 7: The left figure gives m p with correct reweighting factors (blue squares) and with randomly permutedreweighting factors (green diamonds). The right figure is the same but for f p . example parameter points. Note that in both Figs. 7 and 8 we observe an increase in statisticalerrors which appears roughly consistent with what should be expected from the decrease in √ N eff .We should emphasize that further careful studies may be needed to establish a more accurateestimate of possible errors in the reweighting procedure.8 Unitary f π versus m h from reweighting m h total (GeV) f π ( G e V ) Unitary f K versus m h from reweighting m h total (GeV) f K ( G e V ) FIG. 8: Reweighted results for f p (left) and f K (right) as functions of m h at three parameter sets ( b , m l ) :green diamonds: (2.25, 0.008), red circles: (2.13, 0.005), blue squares: (2.25, 0.004).Volume ( m l , m h ) Total MD time Measurement range Measurement total24 (0.005, 0.04) 0-8980 900-8980 every 40 20324 (0.01, 0.04) 1455-8540 1460-8540 every 40 17832 (0.004, 0.03) 0-6756 520-6600 every 20 30532 (0.006, 0.03) 0-7220 1000-7220 every 20 31232 (0.008, 0.03) 0-5930 520-5540 every 20 252TABLE IV: Summary of the five ensembles used in this work. III. UPDATED RESULTS FROM THE ENSEMBLES
In this section we update the results presented on the 24 ensembles in [1] to the extended data setdescribed in Sec. II, and in Table I in particular. For this extended data set we make measurementsof pseudoscalar quantities on a total of 203 configurations for the m l = .
005 ensemble and 178configurations for the m l = .
01 ensemble. These configurations were separated by 40 trajecto-ries as documented in the first two rows of Table IV. In our previous work we used 92 of thesemeasurements on each ensemble [1, 4]. Before performing the analyses we binned the data intoblocks of either 80 or 400 trajectories and the measurements from each bin were then treated asbeing statistically independent. No statistically significant increase in the error was observed withthe analysis using bins of 400 trajectories compared to that with bins of 80 trajectories.9In the following sections the results from the 24 lattices, combined with those obtained on the32 ensembles, will be input into global chiral and continuum fits in order to determine physicalquantities; here we simply tabulate the fitted pseudoscalar masses and decay constants as obtaineddirectly from the correlation functions at our simulated quark masses. In addition, since we usethe mass of the W baryon in the definition of the scaling trajectory, we also present the results for m hhh here together with those for the Sommer scale r and also the scale r . Finally, in Sec.III Awe give the results for the masses of the nucleons and D baryons from the 24 ensembles, althoughthe chiral and scaling behaviour of these masses will not be studied in this paper. We present thesebaryon masses partly for completeness and partly to share our experience in the use of differentsources.On the 24 lattices discussed in this section, the measurements are presented for the two valuesof the sea light-quark mass, m l = .
005 and 0.01, and for the full range of valence quark masses m x , y = . , . , . , . , .
03 and 0 .
04. The ensembles with m l = .
02 and 0.03, presentedin [1], are not included in this paper because such values of m l were found to be too large forSU(2) chiral perturbation theory to describe our data. The value of the sea strange-quark mass inthese simulations is m h = .
04. After completing the global chiral and continuum fits describedin Section V below, we find that the physical value of the bare strange-quark mass, obtained usingthe chiral perturbation theory ansatz, is m s = . ( ) . In this section we anticipate this resultand use reweighting to obtain results also at this value of the strange-quark mass.For the 24 ensembles, we placed Coulomb gauge-fixed wall sources at t = t =
57. Fromeach source, we calculated two quark propagators, one with periodic and the other with anti-periodic boundary conditions. From the periodic propagators for the two sources, denoted by D − P , and D − P , , and the anti–periodic propagators, written as D − A , and D − A , , we form the combinations D − P + A , = (cid:16) D − P , + D − A , (cid:17) and D − P + A , = (cid:16) D − P , + D − A , (cid:17) . (17)The use of periodic plus anti-periodic boundary conditions in the time direction doubles the lengthof the lattice in time, which markedly reduces the contamination from around-the-world propaga-tion in the time direction. For two point functions, such as the propagator of a pseudoscalar mesongiven by h p ( t ) p ( ) i = (cid:229) ~ x Tr (cid:26)h D − P + A , ( t ,~ x ) i † D − P + A , ( t ,~ x ) (cid:27) , (18)on a lattice of time extent N t the time dependence of the contribution of the ground state is given0by h p ( t ) p ( ) i = A [ exp ( − m p ( t − )) + exp ( − m p ( N t − ( t − ))] . (19)Here A is a t -independent constant. For our 24 ensembles, we find that around-the-world propa-gation is not visible in two-point functions. This is not the case however, for three-point functions,as we now explain (although we do not analyze three-point functions in this paper, they are beingevaluated in the computation of B K , for example [34]).For three-point functions of the form h P ( x ) O ( y ) P ( z ) i , where P ( x ) and P ( z ) are pseudoscalar in-terpolating fields and O ( y ) is an operator whose matrix element we wish to measure, we use thewall source at t = P ( z ) and the wall source at t =
57 as the source for P ( x ) . Weonly consider y in the range 5 ≤ y ≤
57, so we do not perform any measurements in the doubledlattice. The doubling of the lattice is important to reject around-the-world propagation in time forsuch measurements. For kaons, we found that a time separation of 52 between the sources gave usa broad plateau, with sufficiently small errors. This measurement strategy was chosen to optimisethe measurement of the kaon bag parameter [4, 34].Before presenting our results for masses, decay constants and r and r , we discuss the values ofthe residual mass and the renormalization constant of the local axial current. The residual mass m ′ res ( m f ) at each partially quenched valence mass used in this work is measured using the ratio [79] m ′ res ( m f ) = h | J a q | p ih | J a | p i , (20)where J a q is the usual DWF mid-point pseudoscalar density composed of fields of each chiralitystraddling the mid-point in the fifth dimension, and J a is the physical pseudoscalar density at thesurfaces of the fifth dimension composed of surface fields in the fifth dimension. The results aregiven in Table V. For completeness we also present the corresponding residual masses obtainedafter reweighting to the physical strange mass in Table VI. The residual mass in the two-flavorchiral limit m res = m ′ res ( m x = m l = ) is given in Table VII and in the left-hand plot of Figure 9.We define Z A to be the renormalization constant of the local axial current, A m , composed of thephysical surface fields. Here we have determined Z A through two methods. In the first, Z A isdetermined for each valence mass using the improved ratio [35] of the matrix element h A ( t ) P ( ) i to h A ( t ) P ( ) i , where A m is the conserved DWF axial current and the results are presented inTable VIII. This method assumes Z A =
1, and we find Z A = . ( ) in the two-flavor chirallimit with the simulated sea strange mass, and Z A = . ( ) when reweighted to the nearbyphysical strange mass. This determination of Z A is illustrated in the plots of Figure 10. As pointed1 m x m l m ′ res ( m x ) measured on the 24 ensembles at the simulated strange quark mass m h = . m x m l m ′ res ( m x ) on the 24 ensembles at the physical strange quark mass. out in [1], we expect Z A = + O ( a m res ) , and in [1] we added a ∼
1% error to account for thesize of this correction. As part of our current work, we have investigated the consequences ofthis correction, which is discussed in detail in appendix B. From this analysis, we find Z A = . ( ) , a 1.8% difference from the result with our previous method. Although, as we will see,this error is smaller than our current combined errors on the decay constants and other physicalquantities, we choose to use this value of Z A = . ( ) , coming from Z V / Z V as defined inEquation (B19), as the normalization factor for the local axial current when quoting all our centralvalues below. Here V and V are the local and conserved vector currents.We now turn to the measurements of the meson masses and decay constants. In order to illustratethe quality of the fits, we start by presenting some sample plots for the unitary pion and kaon onthe m l = . m h = .
04 ensemble. The pion effective masses obtained using different sources2 m h m res m res m sim h m phys h m res in the two-flavor chiral limit on the 24 and 32 ensembles at the simulated and physicalstrange sea-quark masses. m l m r e s ’ m l m r e s ’ FIG. 9: Chiral extrapolation of the unitary values of m ′ res for the 24 (left) and 32 (right) ensembles.While the fit is only marginally acceptable for the 32 lattices, an additional uncertainty of O(5 × − ) isnegligible. and sinks are shown in Figure 11. The mass and decay constant is obtained from a simultaneous fitwith a single, constrained mass to five correlation functions. These are the h P | P i , h A | A i and h A | P i correlation functions (denoted in the figure by PP, AA and AP respectively) with gauge-fixed wallsources and local (LW) or wall (WW) sinks (we do not use the AA-WW combination becauseit is noisier). The long time extent N t =
64 on our lattices together with the noise properties ofpseudoscalar states allow for long plateaux and the results are insensitive to the choice of t min , thestarting point of the fits. Figure 12 displays the effective masses for the unitary kaon, together withthe results obtained from a simultaneous constrained fit. We give an example of the m h dependenceof the unitary pion and kaon masses in figure 13. This dependence is obtained by reweighting.We normalize the states so that, for periodic boundary conditions, the time dependence of thecorrelators for large times is given by C s s O O ( t ) = h | O s | p ih p | O s | i m xy V h e − m xy t ± e − m xy ( N t − t ) i , (21)where the superscripts specify the type of smearing and the subscripts denote the interpolatingoperators. The sign in the square brackets in Eq. (21) is + for PP and AA correlators and − for AP m h Z A (chiral) Z A ( m l = . ) Z A ( m l = . ) m sim h = .
04 0.71651(46) 0.71732(14) 0.71783(15) m phys h Z A on the 24 ensembles at the simulated and physical strange sea-quark masses. t Z A -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 m l Z A FIG. 10: Measurement of Z A for m f = .
005 on the m l = . m h = .
04 ensemble (left panel) and theunitary chiral extrapolation of Z A for the 24 ensembles (right panel). The results do not change significantlyunder reweighting to the physical strange mass. ones. We therefore define the amplitude of the correlator to be N s s O O ≡ h | O s | p ih p | O s | i m xy V . (22)For each correlator included in the simultaneous fit N LWAA , N LWPP , N LWAP , N WWPP and N WWAP , we determine the amplitude and obtain the decay constant f xy using f xy = Z A s m xy N LWAP N WWPP . (23)Table IX contains the measured pseudoscalar masses and decay constants at the simulated strange-quark mass m h = .
04. After reweighting to the estimated physical strange-quark mass m s = . ( ) the masses and decay constants of the pions are presented in Table X and those for thekaons in Table XI.The W baryon, being one of the quantities included in the definition of our scaling trajectory (seeSection V), plays an important rˆole in our analysis. We have performed measurements on the same4 t m xy e ff PP LW t m xy e ff PP WW t m xy e ff AP LW t m xy e ff AP WW t m xy e ff AA LW
FIG. 11: Effective pion masses from the PP LW correlator (top left), PP WW correlator (top right), AP LWcorrelator (center left), AP WW (center right) and AA LW correlator (bottom). Note the different verticalscale for the WW correlators. The horizontal bands represent the result for the mass from a simultaneousfit. configurations using a gauge-fixed box source of size 16 lattice units that gives a good plateaufor the W -state for valence quark masses m x = .
04 and m x = .
03 to enable interpolation to thephysical strange-quark mass. We display the fit to the m x = . W baryon mass on the m l = . m h = .
04 ensemble in figure 14, along with the dependence of this mass on the dynamical strange5 t m xy e ff PP LW t m xy e ff PP WW t m xy e ff AP LW t m xy e ff AP WW t m xy e ff AA LW
FIG. 12: Effective kaon masses from the PP LW correlator (top left), PP WW correlator (top right), AP LWcorrelator (center left), AP WW (center right) and AA LW correlator (bottom). Note the different verticalscale for the WW correlators. The horizontal bands represent the result for the mass from a simultaneousfit. mass using reweighting.The results for the W mass, m hhh , obtained directly at the simulated strange-quark mass ( m h = . m y = .
04 and 0.03 are presented in Table XII. In this table wealso present the results for m hhh obtained after reweighting to the physical strange-quark mass. InTable XIII we display the values of the Sommer scale r , r and their ratio at both the simulatedand physical strange-quark masses. These quantities were determined using Wilson loops formed6 m h m xy m h m xy FIG. 13: We illustrate the m h dependence of the unitary pion (left panel) and kaon (right panel) masseson the m l = . ensemble. The values are obtained by reweighting around the simulated value( m h = . t m xxx m h m xxx FIG. 14: Fit to the W baryon mass with valence strange mass m x = .
04 on the m l = . m h = .
04, 24 ensemble showing the quality of the fit with our box source (left panel). We also show the weak dependenceof the W baryon mass with fixed valence mass m x = .
04 on our simulated m h inferred by the reweightingprocedure on the m l = . ensemble (right panel). from products of temporal gauge links with Coulomb gauge-fixed closures in spatial directions,with an exponential fit to the time-dependence of the Wilson loop W ( r , t ) from t = t = r . The resulting potential V ( r ) was then fit over the range r = . − V ( r ) = V − a r + s r , (24)where V , a and s are constants. These fits are illustrated in Figure 15, which shows the fit to thetime dependence of the Wilson loop W ( r = . , t ) at the physical strange-quark mass, and alsothe subsequent fit over the potential. The strange-quark mass dependence of the scales r and r is7 m x m y m xy ( . ) m xy ( . ) f xy ( . ) f xy ( . ) m xy ( m l ) and decay constants f xy ( m l ) on the 24 ensembles at the simulatedstrange-quark mass ( m h = . small and cannot be resolved within our statistics. A. Nucleon and D Masses
A detailed study of the baryon mass spectrum, including the continuum and chiral extrapolations,is postponed to a separate paper. The one exception is the W baryon, whose mass is used in the8 m x m y m xy ( . ) m xy ( . ) f xy ( . ) f xy ( . ) m xy ( m l ) and decay constants f xy ( m l ) on the 24 ensembles at the physical strange-quark mass m s = . ( ) . m x m xh ( . ) m xh ( . ) f xh ( . ) f xh ( . ) m xh ( m l ) and decay constants f xh ( m l ) on the 24 ensembles at the physical strange-quark mass m s = . ( ) . definition of the scaling trajectory and which is therefore studied in detail together with the prop-erties of pseudoscalar mesons. In this subsection we briefly discuss our experiences in extractingthe masses of the nucleons and D -baryons using different sources and present the results for thesemasses on each ensemble, starting here with those from the 24 ensembles. The baryon spec-trum from the 32 ensembles will be discussed in Sec. IV A. We start however, with some generalcomments about our procedures which are relevant to both sets of ensembles.We use the standard operator, N = e abc ( u Ta C g d b ) u c , to create and annihilate nucleon states and D = e abc ( u Ta C g m u b ) u c for the flavor decuplet D states. On an anti-periodic lattice of size N t in thetime direction, the zero-momentum two-point correlation function, C ( t ) , calculated with one ofthese baryonic operators at its source and sink, takes the following asymptotic form for sufficientlylarge time, t , C ( t ) = Z [( + g ) e − Mt − ( − g ) e − M ( N t − t ) ] , (25)corresponding to particle and anti-particle propagation, respectively. Conventionally one choosesan appropriate range in time where the excited-state contributions can be neglected so that this9 m y m h m W ( . ) m W ( . ) ensembles at the simulated strange quark mass m h = . m h = . m h = . Q ( . ) Q ( . ) Q ( . ) Q ( . ) r r r / r r , r and r / r at the simulated ( m h = .
04) and physical ( m h = . ensembles. Q ( m l ) denotes the quantity measured with light-quark mass m l . form is valid, and extracts the ground-state mass, M , by fitting the numerical data to the functionin Eq. (25). This is indeed what we do to extract baryon masses from the 24 ensembles. Alter-natively we can try to fit the correlation function to a sum of two exponentials, representing theground- and excited-state contributions. As will be reported in Sec. IV A, this is the method weuse for the 32 ensembles.The determination of baryon masses can be made more effective by an appropriate choice ofsmearing at the source and/or sink. We use several different choices of the smearing of theseoperators, wall, box, and gauge-invariant Gaussian [37, 38], in an attempt to obtain a better overlapwith the ground state; our choices are summarized in Table XIV. The wall source, used for the32 ensembles, is Coulomb-gauge fixed. A box source of size 16, also Coulomb-gauge-fixed, isused for the 24 ensembles. The Gaussian-source radius is set to 7 lattice units and 100 smearingsteps are used for the 24 ensembles, while the radius is 6 in the 32 ensembles: these choices areoptimized for our nucleon-structure calculations [7–9].As can also be seen from the table, several steps are taken to reduce the statistical error. For eachconfiguration, as many as four different time slices are used for the sources, usually separated0 t V e ff (r = . ) r V (r) FIG. 15: The effective potential of the Wilson loops with a spatial extent of r = .
45 on the 24 , m l = . t = − V ( r ) on this ensemble, again at the physical strange-quark mass,as a function of the spatial extent of the Wilson loops, overlaid by the fit to the Cornell form over the range r = . − by 16 lattice units, but occasionally fewer. Measurements are made as frequently as every tenthtrajectory and are averaged into bins of 40 hybrid Monte Carlo time units.We now turn to the results obtained specifically on the 24 ensembles. The unitary nucleon and D effective masses are plotted in Figs. 16 and 17 for each choice of quark mass. For the nucleon,both Gaussian and box sources are shown. Plateaus for the effective masses obtained with thebox source appear quickly, suggesting a strong overlap with the ground state. The correspondingplateaus obtained with the Gaussian source appear more slowly, from above. Both sets of resultsagree reasonably well for sufficiently large t . For the D the correlators were only computed usingthe box source and the plateaus for the effective masses again appear quickly. The results forthe masses, obtained using fully correlated fits, are summarized in Table XV. Note such fullycorrelated fits work well for extracting baryon masses as the procedure involves much shorterranges in time than for the meson observables discussed in the rest of this paper. As expected fromthe effective mass plots, nucleon masses obtained using different sources agree fairly well whenthe fits are performed over appropriate ranges. All values of c /d.o.f. are close to 1 or smaller,except for the box-source nucleon fit at m f = .
02 which is about 2.5.Some of these results have been reported earlier at Lattice 2008 [39], and also partially in relatedpapers on nucleon structure [8, 9]. A preliminary report on a bootstrap correlated analysis withfrozen correlation matrix was presented at Lattice 2009 [40] and the results agree with the updated1 size m l source type correlators source time slices configurations24 N D , W N D , W N D , W N D , W N , D
10, 26, 42, 58 2640.004 Wall N , D
0, 16, 32, 48 3050.006 Wall N , D N , D
10, 26, 42, 58 1690.008 Wall N , D
0, 16, 32, 48 254TABLE XIV: Summary of the configurations used in the calculation of the baryon spectrum. m l N (Gaussian) N (Box) D (Box)0.005 0.671(4) { } { } { } { } { } { } { } { } { } { } { } { } TABLE XV: Baryon mass in lattice units from the b = .
13, 24 ensembles. {} denotes fit range. ones given here. IV. RESULTS FROM THE ENSEMBLES
The results for masses, decay constants, r and r obtained directly on the 32 lattice are pre-sented in the same format as those from the 24 ensembles in Section III and the availablemeasurements are also summarised in table IV. The results are presented for three values ofthe sea light-quark mass m l = . , .
006 and 0.008 which correspond to unitary pion masses2 N e ff ec ti v e m a ss Gaussian SourceBox Sourcem l = 0.005 N e ff ec ti v e m a ss Gaussian SourceBox Sourcem l = 0.01 N e ff ec ti v e m a ss Gaussian SourceBox Sourcem l = 0.02 0 2 4 6 8 10 12 14 160.80.840.880.920.961 N e ff ec ti v e m a ss Gaussian SourceBox Sourcem l = 0.03 FIG. 16: Nucleon effective mass plots from the 24 ensembles. Results obtained using the Gaussian sourceare marked by red squares and those from the box source by blue circles. The four plots correspond tounitary light-quark masses 0.005 (top-left), 0.01 (top-right), 0.02 (bottom-left) and 0.03 (bottom-right). in the range 290 MeV – 400 MeV which we had found to be consistent with SU(2) chiral per-turbation theory on the 24 lattice [1]. The valence-quark masses used in the analysis are m x , y = . , . , . , . , .
025 and 0 .
03. For pseudoscalar quantities we use 305, 312and 252 measurements separated by 20 trajectories on the 0.004, 0.006 and 0.008 ensembles re-spectively (see Table IV). For the 32 lattices, we have used a single-source technique for ourmeasurements of pseudoscalar quantities, which differs from the two-source method for the 24 ensembles. Recall that for the 24 ensembles, as discussed in Section III, we placed Coulombgauge-fixed wall sources at t = t =
57. For the 32 ensembles we have used a sin-gle source and calculated both periodic and anti-periodic propagators from this one source. Thesource is placed at t = t src = n mod 64where n is the measurement index, which starts from zero. Moving the source in this way helpsto decorrelate measurements. We always place the anti-periodic boundary condition on the linksin the time direction going from the hyperplane with t = t src − t = t src . Clearly the number ofpropagators to calculate for the single source method is half that for the two-source method.3 D e ff ec ti v e m a ss Box Sourcem l = 0.005 D e ff ec ti v e m a ss Box Sourcem l = 0.010 2 4 6 8 10 12 14 160.920.9611.041.08 D e ff ec ti v e m a ss Box Sourcem l = 0.02 D e ff ec ti v e m a ss Box Sourcem l = 0.03 FIG. 17: Effective mass plots for the D baryon from the 24 ensembles. The results were obtained using thebox source. The four plots correspond to unitary light-quark masses 0.005 (top-left), 0.01 (top-right), 0.02(bottom-left) and 0.03 (bottom-right). For meson two-point functions, as given in Eq. (18), the single-source method is identical to thetwo-source method, except for having half the number of measurements per configuration. For thelight-quark masses on our 32 ensembles we do see around-the-world effects at the fraction of apercent level, so fits of the form in Eq. (19) must be used. We also perform measurements usingthree-point functions of the type h P ( x ) O ( y ) P ( z ) i , where P ( x ) and P ( z ) are pseudoscalar interpo-lating fields and O ( y ) is an operator whose matrix element we wish to measure. Here P ( x ) is madeout of propagators of the form D − P + A , = / (cid:16) D − P , + D − A , (cid:17) in the notation of Eq. (17) and P ( z ) is composed of D − P − A , = / (cid:16) D − P , − D − A , (cid:17) propagators. This means that the time separation be-tween P ( x ) and P ( z ) is N t , the time extent of our lattice. We performed tests on our 24 ensembles,comparing the single-source and two-source methods and found that, for the same number of in-versions, the single-source methods gave at least as small an error as the two-source methods. Thesingle-source method allows us to measure on more configurations for the same computer timeand so we chose this method. Although we do not discuss three-point measurements in this paper,sharing propagators between them and the two-point measurements discussed here has helped todefine our measurement strategy.4 m x m l m ′ res on the 32 ensemble set at the simulated strange quark mass m h = . m x m l m ′ res on the 32 ensemble set at the physical strange quark mass. The measured values of the residual mass m ′ res at each pair of valence and sea light-quark masses( m x , m l ) used in this work are given in table XVI; in this table the strange-quark mass is theone used in the simulation m h = .
03. Table XVII contains the corresponding results obtainedafter reweighting to the physical strange mass ( m s = . ( ) ) determined later in the analysisand presented in Section V. The residual mass in the unitary two-flavor chiral limit is given intable VII and figure 9.The results for Z A for the 32 ensembles obtained from the ratios of matrix elements of A and A are given in table XVIII. We obtain Z A = . ( ) in the chiral limit with the simulated seastrange mass and Z A = . ( ) when reweighted to the nearby physical strange mass. This isillustrated in figure 18. As explained in Section III and appendix B however, in this paper we use Z V / Z V = . ( ) as the normalization factor for the local axial current when calculating the5 m h Z A ( chiral ) Z A ( m l = . ) Z A ( m l = . ) Z A ( m l = . ) m sim h = .
03 0.74475(12) 0.745053(54) 0.745222(45) 0.745328(48) m phys h Z A on the 32 ensembles at the simulated and physical strange sea-quark masses. t Z A m l Z A FIG. 18: Measurement of Z A for m f = .
004 on the m l = . m h = .
03 ensemble (left panel) andthe unitary chiral extrapolation of Z A for the 32 ensemble set (right panel). The results do not changesignificantly under reweighting to the physical strange mass. central values of physical quantities.In order to illustrate the quality of the fits, we present sample effective mass plots for the unitarysimulated pion on the m l = . m h = .
03 ensemble in figure 19 and for the kaon in Figure 20.The analysis is performed as a simultaneous constrained fit to the five pseudoscalar channels asfor the 24 ensembles (see Section III). The fits are performed between t min =
12 and t max = m h dependence of the unitary pion and kaon masses infigure 21.Table XIX contains the measured pseudoscalar masses and decay constants at the simulatedstrange-quark mass m h = .
03. Reweighting to the estimated physical strange-quark mass m h = . ( ) , we obtain the masses and decay constants of the pions and kaons in Tables XXand XXI respectively.We use a gauge fixed box source of size 24 for the W baryon using the same configurations as forour pion measurements with valence strange-quark masses m x = .
03 and m x = .
025 to enablean interpolation to the physical strange-quark mass. We display the fit to the m x = . W baryonmass on the m l = . m h = .
03 ensemble in figure 22, along with the dependence of this mass6
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff PP LW
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff PP WW
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff AP LW
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff AP WW
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff AA LW
FIG. 19: Effective pion masses from the PP LW correlator (top left), PP WW correlator (top right), AP LWcorrelator (center left), AP WW (center right) and AA LW correlator (bottom). Note the different verticalscale for the WW correlators. The horizontal bands represent the result for the mass from a simultaneousfit. on the dynamical strange mass under reweighting. We take our fitting range between t min = t max = W baryon and the scales r , r and r / r are given in Table XXIIand XXIII respectively. r and r were determined again using Wilson loops formed from prod-ucts of temporal gauge links with Coulomb gauge-fixed closures in spatial directions, with an7
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff PP LW
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff PP WW
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff AP LW
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff AP WW
12 16 20 24 28 32 36 40 44 48 52 t m xy e ff AA LW
FIG. 20: Effective kaon masses from the PP LW correlator (top left), PP WW correlator (top right), AP LWcorrelator (center left), AP WW (center right) and AA LW correlator (bottom). Note the different verticalscale for the WW correlators. The horizontal bands represent the result for the mass from a simultaneousfit. exponential fit from t = t = r = . −
10. An example of the fit to the time dependence of the Wilson loops at the physicalstrange-quark mass is given in Figure 23. This figure also shows the fit to the potential. On theseensembles, the strange-quark mass dependence of r and r can be resolved within the statistics,but remains small.8 m h m xy m h m xy FIG. 21: We illustrate the m h dependence of the unitary pion (left panel) and kaon (right panel) masseson the m l = . ensemble. The values are obtained by reweighting around the simulated value( m h = . t m xxx m h m xxx FIG. 22: We display the fit to the W baryon mass with valence strange mass m x = .
03 on the m l = . m h = .
03, 32 ensemble showing the quality of the fit with our box source (left panel). We also show theweak dependence of the W baryon mass with fixed valence mass m x = .
03 on our simulated m h inferred bythe reweighting procedure on the m l = . ensemble (right panel). A. Nucleon and D Masses
Baryon effective masses from the 32 ensembles are plotted in Fig. 24 and 25. The Gaussian-source correlators give good effective-mass signals, while the wall-source correlators are muchnoisier; indeed it is hard to identify a plateau in effective mass signals from the latter. While fornucleons effective mass signals from the wall-source seem to eventually settle at the same valuesas from Gaussian source correlators, for the D baryons a plateau cannot be identified from the wallsource except for the lightest up/down mass. Nevertheless fully correlated fits using two expo-9 m x m y m xy ( . ) m xy ( . ) m xy ( . ) f xy ( . ) f xy ( . ) f xy ( . ) m xy ( m l ) and the decay constants f xy ( m l ) on the 32 ensembles at thesimulated strange-quark mass ( m h = . nentials to represent the contributions of the ground and first-excited states can be performed forboth the nucleon and D , yielding the results summarized in Table XXIV. In addition to this fully-correlated two-exponential fit, we have tried two other fit methods: uncorrelated and bootstrapcorrelated with frozen correlation matrix [40]. While those earlier analysis were conducted onsmaller statistics, they agree with the two-state fully correlated fits within two standard deviations(see Table XXV.) We use the results from the two-state fully correlated fits as our best values of0 m x m y m xy ( . ) m xy ( . ) m xy ( . ) f xy ( . ) f xy ( . ) f xy ( . ) m xy ( m l ) and decay constants f xy ( m l ) computed on the 32 ensembles at thephysical strange-quark mass m h = . ( ) . m x m xh ( . ) m xh ( . ) m xh ( . ) f xh ( . ) f xh ( . ) f xh ( . ) m xh ( m l ) and decay constants f xy ( m l ) on the 32 ensembles at the physicalstrange-quark mass m h = . ( ) . the baryon masses. They also broadly agree with an independent analysis of baryon masses fromour ensembles by the LHP collaboration [41] within two standard deviations. V. COMBINED CONTINUUM AND CHIRAL FITS
We now turn to the main objective of this paper which is to use the results obtained on the 24 and 32 ensembles, as discussed in the previous two sections, to determine physical hadron andquark masses and mesonic decay constants in the continuum limit, for physical values of the lightand strange quark masses. Since we are reporting our first results obtained at a second lattice1 m y m h m W ( . ) m W ( . ) m W ( . ) ensembles at the simulated strange quark mass m h = . m h = . m h = . Q ( . ) Q ( . ) Q ( . ) Q ( . ) Q ( . ) Q ( . ) r r r / r r , r and r / r at the simulated ( m h = .
03) and physical ( m h = . ensembles. Q ( m l ) denotes the quantity measured with light-quark mass m l . spacing, we present a careful discussion of our approach to taking the continuum limit and therelation between evaluating the continuum limit and determining the physical quark masses. Westart in Section V A with a discussion of what we mean by a scaling trajectory and explain insome detail the choice of scaling trajectory which we use in the following. In Section V B wedescribe our power counting scheme, in which we treat the O ( a ) terms in our two ensembles andthe NLO terms in SU(2) chiral perturbation theory as being of comparable size. In order to gaininsights into the uncertainties associated with the chiral extrapolation, in addition to SU(2) chiralperturbation theory, we introduce an analytic ansatz which is a simple first-order Taylor expansionin the light-quark mass. This is explained in Section V C. We then discuss the specific fittingprocedure which implements this power counting strategy in Section V D and in Section V E wepresent and discuss the results. A. Defining the scaling trajectory
Although ultimately we will combine the continuum and chiral extrapolations by performing global fits as described in subsection V A 3 and in the following subsections, we start by focussing2 t V e ff (r = . ) r V (r) FIG. 23: The effective potential of the Wilson loops with a spatial extent of r = .
45 on the m l = . t = − V ( r ) on this ensemble, again at the physical strange-quark mass,as a function of the spatial extent of the Wilson loops, overlaid by the fit to the Cornell form over the range r = . − N E ff ec ti v e M a ss Gaussian SourceWall Source Gaussian 2-exp fitGround-state mass from Gaussian 2-exp fitm l = 0.004 N E ff ec ti v e M a ss Gaussian SourceWall SourceGaussian 2-exp fitGround-state mass from Gaussian 2-exp fitm l = 0.006 N E ff ec ti v e M a ss Gaussian SourceWall SourceGaussian 2-exp fitGround-state mass from Gaussian 2-exp fitm l = 0.008 FIG. 24: Nucleon effective mass plots from the 32 ensembles. on the approach to the continuum limit and discussing the definition and choice of scaling tra-jectory. For the purposes of this subsection we imagine that we can perform lattice computationsfor any choice of quark masses and envision performing a series of lattice simulations for a rangeof values of b , the inverse square of the bare lattice coupling. As b → ¥ the lattice spacing,3 D e ff ec ti v e m a ss Gaussian SourceWall sourceGaussian 2-exp fitGroud-state mass from Gaussian 2-exp fitm l = 0.004 D e ff ec ti v e m a ss Gaussian SourceWall sourceGaussian 2-exp fitGroud-state mass from Gaussian 2-exp fitm l = 0.006 D e ff ec ti v e m a ss Gaussian SourceWall sourceGaussian 2-exp fitGroud-state mass from Gaussian 2-exp fitm l = 0.008 FIG. 25: D effective mass plots from the 32 ensembles. m l N D { } { } { } { } { } { } TABLE XXIV: Nucleon and D masses in lattice units from the 32 ensembles obtained by two-exponentialcorrelated fits to Gaussian-source correlators. {} denotes fit range. measured in physical units, will vanish along with all discretization errors. We refer to such aone-dimensional path through the space of possible lattice theories as a scaling trajectory. For2+1 flavor QCD we must vary the bare lattice mass m ud ( b ) of the up and down quarks and m s ( b ) of the strange quark so that this trajectory describes physically equivalent theories up to order a errors. The functions m ud ( b ) and m s ( b ) can be determined by requiring two mass ratios (or twoother dimensionless quantities) to remain fixed as b varies. Because of the presence of O ( a ) discretization errors, using a different pair of mass ratios will yield a different trajectory of latticetheories, whose low-momentum Green’s functions will be equivalent to those of the first up to O ( a ) corrections.In ref. [1], where we obtained results from simulations at a single value of b , we found that using4 m l full corr. uncorr. bootstrap a LHP b ensembles. Su-perscript a denotes Ref. [40], where a frozen correlation matrix was used and superscript b denotes Ref. [41]. the masses of the p and K mesons and the W baryon to determine the lattice spacing a and the barevalues of m ud and m s was an effective procedure. A natural choice of scaling trajectory wouldtherefore be to keep the ratios m p / m W and m K / m W fixed as b varies. Thus these ratios wouldbe chosen to take their continuum values for all b with no a corrections. This choice of scalingtrajectory then fixes the functions m ud ( b ) and m s ( b ) . In addition, we will identify an inverselattice spacing, expressed in GeV, with each point on this scaling trajectory. To do this we use themass of the W − baryon and define 1 / a = . / m W GeV where 1.672 GeV is the physical massof this baryon and m W is the mass of the W − as measured along our trajectory in lattice units.Having defined the scaling trajectory and determined the lattice spacing at each b by fixing theratios m p / m W , m K / m W and the mass of the W baryon to their physical values, we are in a positionto make predictions for other physical quantities. The results obtained at a particular value of b will differ from the physical ones by terms of O ( a ) . We imagine eliminating these artefacts byextrapolating results obtained at several values of b to the continuum limit. In order to discuss thiscontinuum extrapolation it is convenient to introduce some notation. Let us assume that we haveperformed lattice calculations at a series of N values of b , { b e } ≤ e ≤ N corresponding to pointsalong the scaling trajectory defined above (in the present study N = m e f = m f ( b e ) where f = ud or s . On each of the lattices we computea number of physical quantities, e.g. the kaon leptonic decay constant f e K , and our prediction forthe physical value of f K is the value obtained by extrapolating to the continuum limit.Of course, as already mentioned above, the scaling trajectory and the assigned value of the latticespacing at a particular b are not unique. Had we used three different physical quantities to calibratethe lattice at each b and then used the resulting bare quark masses and lattice spacing to compute m p / m W , m K / m W and the mass of the W baryon, we would find results which differed from the5physical ones by terms of O ( a ) . Although there is a choice of the quantities used to define anddetermine the scaling trajectory and the value of the lattice spacing at each b , for a 2+1 flavortheory the number of conditions is always 3 N , where N is the number of different b values used inthe simulations and the factor 3 corresponds to the fact that at each b there are three parameters,the bare masses m ud and m s and the lattice spacing a .In the above presentation we have tried to provide a pedagogical introduction to the determinationof scaling trajectories and chose to decouple issues related to the extrapolations in the mass of thelight quark (chiral extrapolations) from the discussion. Of course, in practice at present we areunable to perform simulations at physical quark masses, i.e. with masses which give the physicalvalues of m p / m W and m K / m W , and so chiral extrapolations are necessary. It will therefore be usefulin the following to discuss the scaling behavior of a general 2+1 flavor theory in which the massesof the pion and kaon differ from those in Nature. Following the conventions defined elsewhere inthis paper, we will use m l and m h for the quark masses in the DWF lattice action which correspondto the usual ud and s quarks, and e m l and e m h for the corresponding multiplicatively renormalizablebare quark masses e m l = m l + m res and e m h = m h + m res specific to the DWF action. In the nextsubsection we review the origin of the a errors as described by the Symanzik effective theory forDWF and in the following subsection present our treatment of scaling for this more general theory.
1. Symanzik effective theory and a → extrapolation Symanzik’s effective theory provides a powerful framework in which to discuss the approach tothe continuum limit. For any finite value of b we expect the low-momentum Green’s functions inour lattice theory to agree with those in a corresponding effective continuum theory. The effectiveaction for this theory contains not only the usual dimension-3 and 4 terms standard in QCD but alsohigher-dimension operators. If the quark masses and the coefficients of these higher-dimensionoperators are properly chosen then the low-energy Green’s functions of the lattice and effectivetheories will agree through O ( a d − ) provided the effective theory includes all necessary terms ofdimension up to and including d . This implies that the low-energy Green’s functions of the latticetheory and the usual continuum theory will differ by the matrix elements of these dimension-5 andhigher operators which of course are not present in the standard continuum theory.For the domain wall fermion calculation presented here the leading corrections come from opera-tors of dimension 6. While the dimension-5 Pauli term q s mn F mn q is present, its chiral properties6imply that it is generated by chirality violation due to propagation between the left and right do-main walls. This same residual breaking of chiral symmetry gives rise to the residual mass m res ,the coefficient of the dimension-3 mass term which remains when the input quark mass is setequal to zero. The largest value for m res found in our current calculation, m res = . ( ) ,is suppressed from unity by more than two orders of magnitude. Since a similar suppression forthis dimension 5 operator is expected, the combination of chiral symmetry and the small value of a L QCD ∼ . O ( a ) .We require that for our choice of scaling trajectory the matrix elements of these O ( a ) Symanzikterms behave as a , allowing a linear extrapolation in a to give the continuum limit. This impliesthat the coefficients of these operators remain reasonably constant along our trajectory. This istypically achieved by varying only b and quark masses along the trajectory so the only variationin the coefficients of these O ( a ) terms comes from the variations in b which are quite small inpresent scaling studies [80].
2. Scaling and the quark masses
In the present calculation we obtain results using a number of light-quark masses, all of whichare significantly larger than the physical quark masses that were used in the introductory remarksabove to describe a physical scaling trajectory in which m p / m W , m K / m W and m W were fixed attheir physical values. However, we can easily generalize our notion of a scaling trajectory toinclude families of choices for the parameters ( b , e m l , e m h ) for which, in an obvious notation, theratios m ll / m hhh and m lh / m hhh are held fixed. In the language used earlier, we require that the N triplets of parameters ( b e , e m e l , e m e h ) , 1 ≤ e ≤ N , lie on the same scaling trajectory if m ll ( b e , e m e l , e m e h ) m hhh ( b e , e m e l , e m e h ) = m ll ( b e ′ , e m e ′ l , e m e ′ h ) m hhh ( b e ′ , e m e ′ l , e m e ′ h ) (26) m lh ( b e , e m e l , e m e h ) m hhh ( b e , e m e l , e m e h ) = m lh ( b e ′ , e m e ′ l , e m e ′ h ) m hhh ( b e ′ , e m e ′ l , e m e ′ h ) (27)for each pair e and e ′ . The ratio of lattice spacings for such a pair would be defined as a e a e ′ = m hhh ( b e , e m e l , e m e h ) m hhh ( b e ′ , e m e ′ l , e m e ′ h ) . (28)The scaling trajectory determines two functions e m l ( b ) and e m h ( b ) , where these bare masses arenon-trivial functions of b . While a portion of their b dependence should reflect their naive mass7dimension, these quantities also carry a logarithmic dependence on a characteristic of the anoma-lous dimension of the mass operator qq in QCD. Thus, even when expressed as dimensionlessratios, e.g. e m l ( b ) / m W and e m h ( b ) / m W , these parameters will have singular continuum limits (infact, the sign of the anomalous dimension of qq is such that these ratios vanish in the continuumlimit).The mass parameters e m l and e m h are short-distance quantities whose definition is free of infraredsingularities. For example, they could be specified by examining high-momentum, infra-red safeGreen’s functions with no need to compute low-energy masses which are dependent upon thelow-energy, non-perturbative behavior of QCD. While the individual masses e m l ( b ) and e m h ( b ) donot have a continuum limit, both the naive and anomalous scale dependence cancels in their ratio e m l ( b ) / e m h ( b ) , which is well-defined in the continuum limit and agrees with the correspondingratio in conventional renormalization schemes, such as RI/MOM or MS.Let us now assume that we have performed lattice calculations at a series of N values of b , { b e } ≤ e ≤ N , corresponding to points along the scaling trajectory defined above. This will de-termine a series of quark masses e m e f = e m f ( b e ) where f = l or h . It is natural to introduce a seriesof factors which relate the lattice spacings and quark masses between these N ensembles. For con-venience, we identify a primary ensemble , and introduce 3 ( N − ) factors relating each ensemble e to the ensemble as follows: R e a = a a e = m hhh m e hhh (29) Z e f = R e a e m f e m e f for f = l or h . (30)Since the ratio e m l / e m h is well-defined in the continuum limit, the corresponding ratio for eachof these ensembles e m e l / e m e h differs from that limit by a term proportional to ( a e ) . This O ( a ) correction represents the discrepancy between our choice of scaling trajectory with m ll / m lh fixedas we vary b and an alternative choice where instead e m e l / e m e h is held fixed. Since these trajectoriesdiffer at O ( a ) , we expect that e m e l e m e h = lim b → ¥ (cid:18) e m l ( b ) e m h ( b ) (cid:19) (cid:0) + c m ( L QCD a e ) (cid:1) . (31)The term proportional to c m arises from the shifts in m ll and m lh caused by the first-order effects ofdimension-6 terms in the Symanzik effective action. While c m must vanish as e m e l → e m e h , we prefernot to write c m as proportional to the difference e m e l − e m e h because of possible non-analytic terms8in the quark masses (e.g. possible logarithms of m e l ) that may appear in the low-energy matrixelements of these dimension-6 operators. If we divide Eq. (31) evaluated for our primary ensemble by the same equation applied to the ensemble e and Taylor expand in the lattice spacing, weobtain the following useful relation between Z e h and Z e l : Z e h = Z e l (cid:16) + c m L h ( a e ) − ( a ) i(cid:17) (32)implying the 2 ( N − ) Z factors associated with the quark masses actually depend on N quantitiesthrough order a (e.g. we can take the ( N − ) Z e l and c m as the independent quantities). Theconstraints implied by Eq. (32) do not simplify the N = Z h and Z l for the alternative pair of parameters Z l and c m .Equation (32) provides an explicit estimate of how scaling violations revise the standard expecta-tion that all quark masses will scale with a common Z factor as the cut-off is varied. As we will seefrom our simulation results presented below, the terms proportional to c m are small and difficult toresolve from zero given our statistical errors.Since we are now using formulae in which the lattice spacing a e appears alone rather than in aratio, e.g. as a e / a e ′ , it may be useful to explain how we intend this is to be determined. It isnatural to start by considering the physical scaling trajectory discussed in Section V A on which m ll / m hhh = m p / m W and m lh / m hhh = m K / m W . For this physical trajectory, the actual value ofthe Omega mass measured in GeV can be used to define the lattice spacing for any point b e onthat trajectory using a e = m e hhh / ( . ( . ) GeV ) . In our present study, in order to reach thephysical trajectory a chiral extrapolation must be performed from the quark masses used in oursimulation. Ultimately of course, when we present results for dimensionful quantities in physicalunits, it will be necessary to perform the chiral extrapolation and this is the subject of the followingsubsections. For the present discussion of scaling it is sufficient simply to imagine that the latticespacing has been determined in this way and this is the most straightforward way of interpretingthe O (( a e ) ) terms appearing in equations in this subsection. We stress however, that even this isnot strictly necessary. We can consider a scaling trajectory defined by fixed, but unphysical, valuesof m ll / m hhh and m lh / m hhh and define the lattice spacing by assigning an arbitrary value to M hhh ,the mass of the hhh baryon on the trajectory in “physical” units, a e ≡ m e hhh / M hhh . While the valueof a e defined in this way depends, of course, on the choice of M hhh , this arbitrariness is simplyabsorbed by a change in constants such as c m in (31). For the discussion in this subsection it is9sufficient to note that such a definition of the lattice spacing is possible in principle, the numericaldetermination of a e does not actually have to be performed.In the analysis to follow we will examine a family of nearby scaling trajectories in which e m l and e m h vary over limited ranges (specifically, e m l varies up to about 0.013 on our coarser latticeand e m h varies by up to 20% around e m s ). Consider two such trajectories, defined by keeping theratios m ll / m hhh and m lh / m hhh fixed along each trajectory, but taking different values on the twotrajectories. Let m ll / m hhh = r ll and m lh / m hhh = r lh on the first trajectory and m ll / m hhh = r ′ ll and m lh / m hhh = r ′ lh on the second. As b → ¥ , the ratio of bare quark masses on the two trajectorieswill approach a limit up to O ( a ) corrections: e m e f ( r ll , r lh ) e m e f ( r ′ ll , r ′ lh ) = lim b → ¥ e m f ( b ) e m ′ f ( b ) ! (cid:0) + d m , f ( L QCD a e ) (cid:1) , (33)where f = l or h , and e m e l ( r ll , r lh ) and e m e h ( r ll , r lh ) ( e m e l ( r ′ ll , r ′ lh ) and e m e h ( r ′ ll , r ′ lh ) ) are the values of thebare quark masses on ensemble e such that m ll / m hhh = r ll and m lh / m hhh = r lh ( m ll / m hhh = r ′ ll and m lh / m hhh = r ′ lh ). The ratios R a = m hhh ( e m l ( r ll , r lh ) , e m h ( r ll , r lh )) / m e hhh ( e m e l ( r ll , r lh ) , e m e h ( r ll , r lh )) and R ′ a = m hhh ( e m l ( r ′ ll , r ′ lh ) , e m h ( r ′ ll , r ′ lh )) / m e hhh ( e m e l ( r ′ ll , r ′ lh ) , e m e h ( r ′ ll , r ′ lh )) each describe the change inlattice scale as the bare coupling changes from b to b e . In the limit of small bare coupling, thischange of scale can be determined entirely from the short-distance part of the theory and must bethe same for our two trajectories up to order a corrections since these two trajectories differ onlyin the choice of quark masses. Thus we can write R a R ′ a = + d a L (cid:16) ( a e ) − ( a ) (cid:17) (34)where we have explicitly represented the fact that each ratio and hence the ratio of ratios mustapproach unity as a e → a . Both the coefficients d m , f and d a will vanish when the primed andunprimed trajectories that are being compared become identical.Taking the ratio of two versions of Eq. (33), one for b e and the other for our primary ensemble b and using Eq. (34), we obtain an expression for the change in the factors Z f between these twotrajectories: Z e f Z e ′ f = (cid:16) + ( d m , f + d a ) L h ( a ) − ( a e ) i(cid:17) . (35)Since the changes in e m l and e m h between these two trajectories which we wish to compare are small,the resulting coefficients d m , f and d a will also be small and we will neglect the O ( a ) correctionon the right-hand side of Eq. (35). Thus, we will use the same values for Z l and Z h for this family0of nearby trajectories, i.e. we drop lattice artefacts proportional to ˜ m l and ( ˜ m h − ˜ m s ) and so neglectthe mass dependence of Z l and Z h in this limited range of masses. In the following we will referto this range for e m l and e m h as their “allowed range”.
3. Fitting strategies
We exploit the above relations between numerical results obtained at the two values of b for whichwe have performed simulations in two ways. The first we label the “fixed-trajectory” method. Inthis approach we determine R a , Z l and Z h by matching results obtained at a single pair of equivalentquark masses [81]. For example, the masses used at one value of b may correspond to values atwhich a simulation was actually performed. The corresponding set of masses for the other b mightbe determined by linear interpolation to make the two ratios m ll / m hhh and m lh / m hhh agree withthose on the first ensemble. The ratio of lattice spacings and the two Z f factors are then determinedfrom Eqs. (29) and (30). It will be important to recall that Z l and Z h are constant in the allowedrange of quark masses. Finally, knowing the three factors R a , Z l and Z h we make a common fit tothe mass dependence of physical quantities computed for both values of b .In the final step, we adopt an ansatz for the mass dependence that is expected to be accurateboth for the points in our calculation and for the physical values to which we wish to extrapolate,specifically a NLO chiral expansion about the chiral limit or a simple Taylor expansion aboutthe physical point. Each ansatz for the continuum theory, when combined with the three scalingfactors R a , Z l and Z h and with any required a corrections, will then provide a set of formulaewhich should describe all of our data for both b values. For example, in the chiral fits described inthe next section we can use a common set of Low Energy Constants (LECs) to fit both sets of dataprovided we scale the values used on one set by the required factors of R a , Z l and Z h before weuse them on the other. Where explicit O ( a ) terms are required, these can be added with unknowncoefficients which are also scaled appropriately between our two values of b . In such a combinedchiral and a expansion we adopt a power counting scheme, described below, so that only effectsof a similar minimum size are consistently included.During the initial process of determining R a , Z l and Z h we cannot assign a physical value to thelattice spacing. The original trajectory being used does not correspond to physical masses so nonotion of “GeV” exists for that case. Of course, the further fitting to the quark mass dependence ofthe two ensembles is introduced to allow extrapolation to physical values for the ratios m ll / m hhh m lh / m hhh . When m W is evaluated at this same physical point, its value can be compared with1.672 GeV to determine the lattice scale.This fixed trajectory method is intended to cover a wider range of possible scaling trajectories thanthe example discussed above where the trajectory passes precisely through one of the simulationpoints. If we wish, we can adopt an ansatz for the quark mass dependence of m p , m K and m W andperform this fixed trajectory scaling with the parameters R a , Z l and Z h allowed to vary and fix theirvalues from Eqs. (29) and (30) at values of m l and m h for which the ratios m ll / m hhh and m lh / m hhh take their physical values.The second approach, termed “generic scaling”, introduces the factors R a , Z l and Z h as parametersinto the ansatz being used to fit the quark mass dependence. In this approach we perform a fit to allour data for m p , m K and m W over a range of quark masses for which the fitting ansatz is accurateand for which the use of fixed values for R a , Z l and Z h is legitimate. In this generic scaling ap-proach, our choice of scaling trajectory with fixed hadron mass ratios m ll / m hhh and m lh / m hhh andwith m hhh determining the lattice scale is realized somewhat indirectly. The three conditions asso-ciated with this choice of scaling trajectory are realized by omitting possible a corrections fromthe expressions used to fit m ll , m lh and m hhh . The resulting trajectory can therefore be interpretedas being the one along which the masses of the pion, kaon and W -baryon take their physical values,as was the case in the discussion of Section V A. The difference of course, is that whereas in Sec-tion V A we envisaged (unrealistically at present) being able to simulate directly at the physicalvalue of m l , we now reach the physical point after an extrapolation in quark masses. The detaileddiscussion of the ChPT functions used in describing the quark mass dependence of the pion andkaon masses is given in Subsection V B and those for the analytic ansatz in Subsection V C below.However, both our ChPT and Taylor expansion ans¨atze stipulate that to the order being studied m hhh is a linear function of e m l and e m h . It is instructive to explore this case here.Included among the equations used to determine the low energy constants and the scaling factors R a , Z l and Z h are two equations for m hhh on our two ensembles: m hhh ( e m l , e m h ) = m hhh ( , e m h ) + c m W m l e m l + c m W m h ( e m h − e m h ) (36) m hhh ( e m l , e m h ) = R a m hhh ( R a Z l e m l , R a Z h e m h )= R a h m hhh ( , e m h ) + c m W m l R a Z l e m l + c m W m h (cid:16) R a Z h e m h − e m h (cid:17)i . (37)Here is our primary ensemble, for us that is the one with b = .
25 and the 32 ×
64 volume, whilethe second ensemble is the one with the coarser lattice spacing and is labeled . m e hhh ( e m l , e m h ) are2the hhh -baryon masses corresponding to bare-quark masses e m l and e m h on ensemble e . Althoughwe have written e m h as a general constant, we have in mind to use the equations with e m h in theallowed range of the physical bare strange quark mass in the primary ensemble. Equations (36)and (37) define the three constants m hhh ( , e m h ) , c m W m l and c m W m h which are related to the physical W − mass and its “physical” dependence on the quark masses. The absence of O ( a ) correctionsto Eqs. (36) and (37) implements our choice that m W is being used to set the scale and hence byconstruction contains no finite lattice-spacing errors. While part of a larger set of equations whichare being used to determine the low energy constants as well as R a , Z l and Z h , the leading ordereffect of these two equations is to determine R a . Note that this is identical to imposing Eq. (29)in the fixed trajectory method at the point e m l = e m h = e m h . Since the variation of R a as e m l and e m h change over their allowed range is of the same size as the variation of Z l and Z h over this samerange it can also be neglected, so any particular choice of e m h is equivalent to any other within thisallowed range.The fixed trajectory and generic scaling methods are similar in nature. Both require that an ansatzbe adopted to allow the quark mass dependence of lattice quantities to be described in order todefine the scaling parameters R a , Z l and Z h and to extrapolate to the physical point. Both assumethat the scaling relations between the two ensembles defined by R a , Z l and Z h hold over the allowedrange of masses. The fixed trajectory method corresponds most closely to our original definitionof a scaling trajectory and decouples the matching of the two lattices from the chiral extrapolation.It requires however, the introduction of a convenient but arbitrary point at which the matchingbetween the two ensembles is performed. The generic method avoids this arbitrary choice andapplies these assumptions uniformly over the entire range of allowed masses. The fixed trajectorymethod determines R a , Z l and Z h in an iterative fashion as explained in Section V D. The genericapproach determines the coefficients in the adopted ansatz from a single c minimization. Thephysical quark masses are then determined by inverting the resulting equations which give m p , m K and m W in terms of e m l and e m h .The detailed discussion and results presented in this paper correspond to the fixed trajectorymethod; fits using the generic scaling approach were performed to monitor the consistency ofthe results and estimated errors.3 B. Scaling and chiral perturbation theory
At the start of section V A we discussed the continuum extrapolation in an idealized situation inwhich we can perform simulations at any value of the quark mass m l . In reality this is not the case;for example, the lightest unitary pion appearing in the current study has mass 290 MeV. In order tocompare our results with Nature we therefore need to extrapolate to lighter quark masses and thiswas already acknowledged when discussing the fitting strategies in section V A 3 above. We nowexplain how we combine the continuum and chiral extrapolations in global fits . We start in thissection by using SU(2) chiral perturbation theory for the mass dependence, with the expectationthat the extrapolation will be made more precise if constrained by the theoretically known behaviorof QCD in the chiral limit [1]. However, in order to estimate possible systematic errors associatedwith this extrapolation and to obtain a more complete understanding of the implications of ourcalculation, we also examine a simpler analytic extrapolation to physical quark masses [42] andthis is explained in the following subsection. Although later we will perform extrapolations usingpartially quenched ensembles, for the purposes of this introduction we restrict the discussion tothe unitary theory in which the valence and sea quark masses are equal.We now explain the power counting scheme we employ to identify NLO corrections to the chiraland continuum limits. Since the pion mass and decay constant are central to SU(2) ChPT, webegin by considering the predictions of continuum NLO ChPT for these two quantities: m ll = c l + c l · ( f (cid:16) ( L ( ) − L ( ) ) + ( L ( ) − L ( ) ) (cid:17) c l + p f c l log c l L c ) (38) f ll = f + f · ( f ( L ( ) + L ( ) ) c l − c l p f log c l L c ) . (39)Here m ll and f ll are the mass and decay constant of the pseudoscalar meson composed of twolight quarks, f , L , L , L and L are the conventional low energy constants and L c is the usualchiral scale. The quantity c l comes directly from the lowest order chiral symmetry breaking termin the effective chiral theory and is proportional to the QCD light quark mass. It is conventionallywritten c l = B e m l , where B is another low-energy constant.We now discuss how we apply these formulae to describe the low energy behavior of lattice the-ories which lie on a scaling trajectory. For a sequence of ensembles { e } ≤ e ≤ N lying on such ascaling trajectory not only will the quark masses and lattice units, ( e m e l , e m e h , a e ) be related, but also,when expressed in physical units, the quantities f , L , L , L and L should take the same val-4ues up to O ( a ) corrections. The same is true for the renormalization independent combination c l = B e m l (see the discussion below). As detailed in Ref. [1], chiral perturbation theory at finitelattice spacing for domain wall fermions involves a simultaneous expansion in the explicit barequark mass, m l , the squared lattice spacing, a , and the residual chiral symmetry breaking arisingfrom the finite separation, L s , between the two four-dimensional walls in the fifth dimension. Wewill denote this last quantity by e − l L s , suggesting the exponential decrease in such residual chiralsymmetry breaking found in perturbation theory for DWF. (The actual behavior is a sum of ex-ponential and inverse power dependence on L s .) No new terms need to be added to the resultingeffective low energy theory to describe the resulting Green’s functions to NLO in the parameters e m l , a and e − l L s . Thus, we can use equations with the form of Eqs. (38) and (39) to describethe lattice results for m ll and f ll along a scaling trajectory, provided we work to NLO in a powercounting scheme which treats the quantities c l / ( p f ) , a L and e − l L s as equivalent and keepa single power of any of these quantities as a correction. We must now determine how the param-eters appearing in these equations must be adjusted to describe lattice results at finite a .Since the scale L c can be freely varied if the other analytic terms are appropriately changed, wewill choose this quantity to be constant if measured in physical units. Thus, for each point onour physical scaling trajectory we will choose L c = m W · / . c l all of the parameters which appear in thelarge curly brackets on the right hand side of Eqs. (38) and (39) can be given their continuumvalues, dropping possible O ( a ) terms as being of NNLO in our power counting scheme. Thus,within those brackets the quantities f , L , L , L and L , when expressed in physical units, can begiven identical values for the ensembles on the scaling trajectory.In contrast, when Eq. (39) is used to describe our finite lattice spacing results, the LO quantity f e determined on ensemble e , expressed in physical units, depends on b e . However, it approaches itscontinuum limit with O ( a ) corrections and so we write f e = f + c f ( a e ) .Given the definition of a scaling trajectory, the variation of the quantity c e l needed to apply Eq. (38)to the ensemble e is actually trivial. Because our choice of quark mass e m e l gives the same value for m ll for each ensemble e on our scaling trajectory, all of the quantities in Eq. (38) with the possibleexception of the c e l which we are now considering, are the same when expressed in physical unitsfor all points on the scaling trajectory. Thus, c e l = B e e m e l / ( a e ) must be a constant as well, where B e and e m e l are explicitly left in lattice units. Since we know how the quantities e m l and a are relatedbetween an ensemble e and our primary ensemble , we can determine the N − B e in5terms of the single constant B : B e = Z e l R e a B (40)without any a corrections. Because of the complex scaling behavior of the mass, we will treat B as one of the LEC’s to be determined in our fitting and not relate it to a “physical” continuumquantity whose definition would require introducing a continuum mass renormalization scheme.We conclude that our lattice results for light pseudoscalar masses and decay constants obtainedfrom a series of ensembles { e } can be described through NLO by the formulae: ( m e ll ) = c e l + c e l · ( f (cid:16) ( L ( ) − L ( ) ) + ( L ( ) − L ( ) ) (cid:17) c e l + p f c e l log c e l L c ) (41) f e ll = f (cid:2) + c f ( a e ) (cid:3) + f · ( f ( L ( ) + L ( ) ) c e l − c e l p f log c e l L c ) (42)with c e l = Z e l R e a B e m e l ( a e ) (43)where all quantities in Eqs. (41) and (42) are expressed in physical units (except for B and e m e l inEq. (43) which are given in lattice units).Two important refinements should be mentioned. First, for the case of a physical scaling trajectory,i.e. one which terminates in the physical masses m p , m K and m W , these physical units are naturallyGeV. However, for other scaling trajectories appropriate “physical” units to use can be those inwhich the Omega mass is unity. Second, for simplicity in Eqs. (38), (39), (41) and (42) we havetreated the heavy quark mass as fixed and not displayed the dependence of the quantities f , B , L , L , L and L on m h . In practice we can easily generalize these equations to describe thedependence of m ll and f ll on m h as well. Provided we limit the variation of m h to a small rangeabout an expansion point e m h , this variation can be described by including a linear term in m h − e m h and treating this term as NLO in our power counting scheme. Thus, such extra linear terms willonly be introduced into the leading order terms in Eqs. (41) and (42).Next we present the corresponding formulae for the quantities m K and m W which are used in thedetermination of the scaling trajectory and in the assignment of a lattice spacing at each value of b : ( m e lh ) = (cid:16) m ( K ) (cid:17) + (cid:16) m ( K ) (cid:17) (cid:26) l + l f c e l (cid:27) (44) m e hhh = m ( W ) + m ( W ) c m W , m l c e l . (45)6Here m ( K ) and m ( W ) are the mass of the lh meson and the hhh baryon respectively in the SU(2)chiral limit, i.e. with e m l =
0, for the value of e m h used in the simulation. Similarly the LECs l , and c m W , m l depend on e m h and we are using the notation for the LECs l , which we introducedin [1]. (Note that c m W , m l , whose value is given in Table XXVII below, should be distinguishedfrom the related parameter c m W m l which appears in Equations (36) and (37) above.) The absenceof any corrections of O ( a ) on the right-hand sides of Eqs. (44) and (45) follows from the sameargument which justified omitting an O ( a ) correction from the right hand side of Eq. (41). Formasses e m e l and e m e h lying on a scaling trajectory the left hand sides of these equations must all bethe same because of our definition of scaling trajectory. Because of our power counting scheme,no a corrections need to be included in the NLO terms proportional to c e l on the right hand sideof these two equations. Therefore the leading order terms m ( K ) and m ( W ) must also be the same forall ensembles when expressed in physical units and no O ( a ) correction can appear. As discussedabove, these equations can be generalized to describe the NLO dependence on e m h varying aboutan expansion point e m h . In fact, for the W baryon this more general case for Eq. (45) was describedin the previous subsection in the equivalent Eqs. (36) and (37).Note that the coefficient of the chiral logarithm in Eq. (41) includes a factor which depends on f ,the pion decay constant in the SU(2) chiral limit (all other factors of f in Eqs.(41) and (44) can beabsorbed into a redefinition of LECs which in any case are determined by fitting). This low energyconstant f can be determined from the measured values of f ll using Eq. (42), but to NLO it canalso be replaced by the measured values of f ll .As described in Subsection V A 3, these ChPT formulae can now be used to determine physicalresults in the continuum limit from those obtained on our two lattice spacings. We can employ thefixed trajectory method, finding the ratios Z l and Z h which relate a specific choice of quark masseson one ensemble to those on the other which lie on the same scaling trajectory. The correspondingratio of values of m hhh determines R a . These three quantities then allow a single set of LECs tobe used to extrapolate the results of both ensembles to the continuum limit and to the physicalvalue of the light quark mass using Eqs. (41), (42), (44) and (45). As a result we learn the physicalvalues of e m ud ( b e ) , e m s ( b e ) and a e on our two ensembles. In other words, we determine the quarkmasses and lattice spacings for our two ensembles which lie on the physical scaling trajectory.Alternatively, we can use the generic fitting approach and introduce the three parameters ( Z l , Z h , R a ) into the four equations Eqs. (41), (42), (44) and (45) and obtain a fit to the lattice datafrom both ensembles for which the quark masses lie in the allowed range. The resulting values of7the LECs and ( Z l , Z h , R a ) then determine the functions m e ll ( e m l , e m h ) , m e lh ( e m l , e m h ) and m e hhh ( e m l , e m h ) .The physical quark masses on each ensemble, m e ud = m ud ( b e ) and m e s = m s ( b e ) , are then obtainedby solving the equations: m e ll ( e m e ud , e m e s ) m e hhh ( e m e ud , e m e s ) = m p m W and m e lh ( e m e ud , e m e s ) m e hhh ( e m e ud , e m e s ) = m K m W , (46)where on the right-hand sides the ratios take their physical values.Having determined m ud ( b e ) , m s ( b e ) and a e as described above, we are in a position to computeother physical quantities. For example, at NLO in our power counting the behaviour of the kaondecay constant f K is f e lh = f ( K ) h + c f ( K ) ( a e ) i + f ( K ) ( l + l f c e l − ( p f ) c e l log c e l L c ) , (47)where f ( K ) is the result in the SU(2) chiral limit ( e m l = l , are m h -dependent low-energyconstants and c f ( K ) is a constant. For each b e , having determined e m s ( b e ) we measure f e lh for e m e h = e m s ( b e ) as a function of e m l ; fit the measured values at all b e to determine the LECs and c f ( K ) in Eq. (47) and finally obtain the physical value of f K by setting a = e m l = e m ud . Such aprocedure is then generalized to the other physical quantities we wish to compute. C. Scaling combined with an analytic ansatz for the chiral dependence
While we know that the ansatz based on chiral perturbation theory described in the previous sub-section is valid in the limit of small u and d quark masses, we do not know the precision withwhich it holds over the range of masses which we analyze in this paper (corresponding to data inthe range 240 MeV ≤ m p .
420 MeV). Indeed it is precise lattice simulations which will answersuch questions. In order to obtain some understanding of the corresponding systematic uncertain-ties, in addition to the procedures based on chiral perturbation theory described in section V B, weconsider an ansatz based on a first-order Taylor expansion about a non-zero quark mass, in thestyle of ref. [42, 43]. Within this approach, since we do not include chiral logarithms, we are notable to take the chiral limit and only assume the validity of the analytic ansatz between the physicalpoint (to which we extrapolate) and the region where we have data. In this work we only considerlinear, first-order fits and are therefore insensitive to the choice of expansion point which we taketo be the same as that at which we match the ensembles when using the fixed trajectory method.This simplifies the discussion below of the simultaneous expansion in a and mass differences.8Beyond first order, convergence may be improved by considering an expansion point between theregion in which we have data and the physical point, but this is beyond the scope of our currentanalysis.Using the analytic ansatz for m p as a function of the quark mass m q , we find numerically that theconstant (mass independent) term is consistent with zero, indicating that the tangent of m p ( m q ) in the unitary case does pass through the origin. Thus, at our statistical precision, no significantchiral curvature is needed to satisfy Goldstone’s theorem, however we retain the view that we areindeed using a model which is valid only in a restricted region of non-zero quark masses.Goldstone’s theorem also applies in the partially quenched theory and the pion mass vanishes asthe valence-quark masses are taken to zero while keeping the sea-quark masses fixed. In this casehowever, our linear fit extrapolates to a non-zero pion mass for massless valence quarks, and thisnaturally implies that some form of curvature is required at smaller masses. This is consistent withenhanced chiral logarithms in the partially quenched theory. However, the fits do not necessarilyimply that chiral logarithms at NLO correctly represent the quark-mass dependence between thesimulated range of masses and the physical point. Instead, in this approach the sum over multipleorders of chiral perturbation theory is assumed to be approximated by a linear dependence in therelevant range of masses. It is also possible of course that the simulated range of masses is outsidethe useful domain of chiral perturbation theory and that, for example, phenomenological modelsbased on combining NLO chiral perturbation theory with arbitrary analytic subsets of terms whichappear at NNLO and NNNLO are less well motivated than our linear ansatz.For m p and f p it is convenient to define the average valence quark mass e m v = e m x + e m y . As insection V B, we apply a power counting rule in a double expansion in m x − m m , m y − m m , m l − m m and a , where m m is the mass at which we match the ensembles which we also choose to be thepoint around which we perform the Taylor expansion and we recall that m x , y and m l are the valenceand sea light-quark masses respectively (here we allow for partial quenching). For the pion masswe use the ansatz m xy = C m p + C m p ( e m v − e m m ) + C m p ( e m l − e m m ) , (48)where we use our standard notation in which the subscripts xy imply that the two valence quarkshave mass m x and m y respectively. By the definition of our scaling trajectory, there is no O ( a ) term at the match point and so there is no correction to C m p . Within our power counting we couldequivalently use m xy = C m p + C m p e m v + C m p e m l , (49)9where for convenience we redefine C m p between equations (48) and (49).In searching for evidence of chiral logarithms it is conventional to plot the ratio m xy / e m v as afunction of the quark masses. With the ansatz proposed in Eq. (49) m xy e m v = C m p e m v + C m p + C m p e m l e m v , (50)and we note that an observed deviation of the mass dependence of m xy e m v from a constant in thefinite range of quark masses which can be simulated, is not in itself unambiguous evidence of anon-analytic structure.For decay constants, which do not vanish in the chiral limit, the O ( a ) term are not sensitive to thechoice of expansion point: f xy = C f p [ + C f p a ] + C f p ( e m v − e m m ) + C f p ( e m l − e m m ) (51) ≡ C f p [ + C f a ] + C f p e m v + C f p e m l , (52)where again we have redefined C f p between the first and second lines.Following a similar argument, at a fixed strange-quark mass, we take the light-quark mass depen-dence of the kaon mass and decay constant and the mass of the W -baryon to be given by m xh ( a , m l ) = C m K + C m K e m x + C m K e m l , (53) f xh ( a , m l ) = C f K [ + C f K a ] + C f K e m x + C f K e m l . (54) m hhh ( a , m l ) = C m W + C m W e m l . (55)We stress that the constants C m p n , C f p n , C f , C m K n , C f K n , C f K and C m W n implicitly depend on the strangequark mass. D. Procedure for combined scaling and chiral fitting
Having introduced the theoretical framework behind our combined scaling and chiral fits in Sec-tions V B and V C we now explain its practical implementation. The formulae given above whichdescribe the combined behaviour are valid only for a fixed strange-quark mass and we are facedwith the problem that the physical strange mass is not known a priori but is an output of thecalculation. The procedure for performing the combined chiral-continuum fits is therefore neces-sarily iterative. As explained in more detail below, we start with some initial values for the lattice0spacings and quark masses, perform the fits and then use linear interpolations in m h to obtainupdated estimates. The process terminates when the updated estimates converge. During this it-erative procedure we use reweighting (see section II D) to adjust all pionic observables to the newstrange-quark mass on each ensemble. For kaon and W observables a linear interpolation betweenthe unreweighted unitary measurement, and measurements with a second valence strange quark(reweighted-to-be-unitary) suffice to obtain that observable for m y = m h = m guess s .For the remainder of this subsection we explain further the procedure which we use to matchlattices with different b and present results for the ratios R e a and Z e f defined in Eqs. (29) and (30)for our ensembles using the fixed trajectory method explained in Section V A 3. We start by takinga specific value of ( m l , m h ) M on the ensemble M to which the other ensembles are matched. Werefer to this as the matching point. The ensemble set M may be the same as the primary ensemble , but does not need to be. As discussed in section V A, the matching to other ensembles e = M isperformed by requiring that the ratios of hadronic masses m ll m hhh and m lh m hhh are the same on all latticesat the matching point. Although the final physical predictions do not depend upon the choice ofmatching point, certain choices are favoured due to the quality of the data at the matching pointand the range over which the data must be interpolated/extrapolated on the other ensembles toperform the matching. The ideal point has as small a statistical error as possible and lies withinthe range of simulated data on all of the matched ensembles such that only a small interpolation isrequired. In practice, the errors on the mass ratios at the matching point can be reduced by fittingto all partially quenched simulated data on the ensemble set M and interpolating to the matchingpoint along the unitary curve. We use linear fitting functions for the light-quark mass dependenceof the pseudoscalar mesons and the W baryon in these short interpolations: m xy = c + c l m l + c v ( m x + m y ) , (56) m xh = d + d l m l + d v m x , (57) m hhh = e + e l m l , (58)where as elsewhere x , y ( l ) represent the light valence (sea) quarks and h represents the heavyquark. Equations (56) - (58) are written in lattice units. Although the linear behaviour in Eqs. (56) -(58) is similar to that used in the analytic ansatz, Eqs. (49), (53) and (55), we stress that themeaning is different. When using the analytic ansatz we assume its validity in the full rangeof masses between the physical ones and those we simulate. Eqs. (56) - (58) on the other hand,are only assumed to represent the mass behaviour in the short intervals between the matching1and simulated points on ensembles e = M , independently of whether we subsequently use chiralperturbation theory or the analytic ansatz to perform the chiral extrapolation.Once a matching point has been chosen, the matching proceeds as follows:1. For each set of ensembles e = M , we perform an independent partially-quenched linear fitto the simulated pion, kaon and Omega masses using the forms given in Eqs. (56) - (58).2. We make a first estimate of the pair of quark masses ( m l , m h ) e on each ensemble set e = M that corresponds to the matching point.3. We then interpolate the three hadronic masses to the estimated m e l for each value of thesimulated unitary heavy quark mass.4. We linearly interpolate each quantity to the estimated value of m e h .5. Next we calculate the ratios R e l = m e ll m e hhh and R e h = m e lh m e hhh .6. Using the measured slopes of m e ll and m e hhh with respect to m e l , by comparing R e l to thecorresponding value R M l at the matching point we obtain an updated estimate of m e l .7. Similarly, by comparing the ratio R e h to R M h we obtain an updated estimate of m e h .8. With these updated estimates of the quark masses ( m l , m h ) e , we return to step 3 and iteratethe steps until the process converges.Once this procedure has converged, we have a set of bare quark masses ( m l , m h ) e which, in phys-ical units, are equivalent to the masses ( m l , m h ) M . Following the discussion in Sec. V A 2, wechoose a primary ensemble and determine the ratios of quark masses Z e f in ensembles and e asin Eq. (30) with the corresponding ratios of lattice spacing R a given in Eq. (29).In the above we assumed that for each ensemble e we had performed simulations at several val-ues of m e h . In our present study the simulations were performed at a single value of m e h and thedependence on the heavy-quark mass is obtained by reweighting as explained in Section II D.The above discussion was deliberately presented in a general case where there are an arbitrarynumber of ensembles. In our case we only have two sets, i.e. the 24 and 32 lattices. For theprimary ensemble we choose the finer 32 lattice. As we have only one other ensemble set (24 ),from now on we drop the superscript on the ratios of lattice spacings ( R a ) and quark masses ( Z l and Z h ).2 M ( am l ) M ( am h ) M ( am l ) e ( am h ) e Z l Z h R a Z l and Z h and the lattice spacing ratio R a determined bymatching at five points over both ensemble sets. The quark masses here are quoted without the additive m res correction. The ensemble e = M . Q / Q r f ll r f lh r m ll r m lh r m hhh f ll /m hhh f lh /m hhh m ll /m hhh m lh /m hhh f ll /m lh f lh /m lh m ll /m lh f ll /f lh m ll /f lh FIG. 26: Ratios of dimensionless combinations of lattice quantities Q (listed in the figure) between the 32 and 24 lattices at the matching point corresponding to m l = . m h = .
03 on the 32 lattice. A value ofunity indicates perfect scaling. The ratios m ll / m hhh and m lh / m hhh (and consequently m ll / m lh ) are definedto scale perfectly at these quark masses as a consequence of our choice of scaling trajectory. In Table XXVI we give results for Z l , Z h and R a obtained by matching at several matching pointson both ensemble sets M ∈ { , } . Since we prefer to have a matching point within the rangeof simulated data on both ensembles, we can discard the first and last entries in the table. Fromthe remaining 3 possibilities, we choose as our final values Z l = . ( ) , Z h = . ( ) and R a = . ( ) from the second entry with M = and ( m l , m h ) = ( . , . ) .Having chosen to perform the matching of the lattices at the two lattice spacings by requiringthat m ll / m hhh and m lh / m hhh take the same values at the matching point, we expect to see latticeartefacts in ratios of other physical quantities. This is illustrated in Figure 26 in which we show theratios of several other dimensionless combinations of lattice quantities between the two lattices at3the quark masses used in the matching procedure above. The figure shows that we can expect onlysmall scaling violations on the order of 1–2% for the other quantities used in our global fits, andalso confirms that other dimensionless combinations of lattice quantities would be equally suitablechoices for the definition of the scaling trajectory. E. Results of combined scaling and chiral fits
Using the matching factors Z l , Z h and R a determined as described in the previous section we areready to perform a simultaneous fit of all our pion, kaon and W mass and decay constant datato either the NLO forms in chiral perturbation theory, Eq. (41) to Eq. (45), or the analytic formsEq. (49) to Eq. (55). We also correct for finite volume effects in NLO PQChPT by substituting thechiral logarithms with the corresponding finite-volume sum of Bessel functions [44]. The iterativeprocedure is the same for each of these three fit ans¨atze. For each iteration i , we:1 estimate the physical strange-quark masses, m is , from the ( i − ) th iteration;2 interpolate and reweight the data to m is ;3 fit the m x , m y , m l dependence of the light pseudoscalar mass and decay constant;4 fit the m x , m l dependence of kaon quantities at m h = m is ;5 fit the m l dependence of the Omega mass for m h = m is ;6 by comparing to the physical values of m p / m W and m K / m W , determine the iterated predic-tions for the physical strange quark masses m i + s .This process is repeated until it converges and a self consistent set of quark masses, lattice spacingsand results in the continuum limit are obtained.For the fits based on NLO chiral perturbation theory we use Eqs. (41) and (42) for the pion massand decay constant respectively, and Eqs. (44) and (47) for the kaon mass and decay constant.In our earlier work [1] we found that we had to apply cuts to keep the pion mass below around420 MeV in order for NLO SU(2) ChPT to give an acceptable description of our data. All theadditional data introduced in this work satisfies this cut and we include all the data for pions withvalence masses m x , m y ≤ .
01 on the two 24 ensembles and all data for pions with valence masses m x , m y ≤ .
008 for the three 32 ensembles. For kaons we include all the valence light-quark4 Parameter No FV Corrections With FV Corrections B f c f L ( ) -0.00000(7) -0.00005(7) L ( ) L ( ) -0.00003(4) -0.00005(4) L ( ) m ( K ) f ( K ) c f ( K ) l l l -0.0018(9) -0.0016(10) l m ( W ) c m W , m l − − TABLE XXVII: Parameters of the global fit to our ensembles using NLO ChPT without finite-volumecorrections (second column) and with finite-volume corrections (third column). For the unitary theory theparameters are defined in Sect. V B and for the partially quenched theory in appendix B of Ref. [1]. masses in the above range for each fixed strange-quark mass. For this infinite-volume SU(2) NLOglobal fit the fitted parameters are presented in the second column of table XXVII. The c / doffor all the fits discussed here are given in table XXVIII. We also perform the corresponding fitsusing the finite-volume chiral logarithm composed of a sum of Bessel functions [44]; resummedexpressions are not available for our partially quenched fits. The parameters of the fit are presentedin the third column of table XXVII. In terms of the conventional LECs ¯ l and ¯ l the results are¯ l = . ( ) , ¯ l = . ( ) ( Infinite Volume ChPT ) (59)¯ l = . ( ) , ¯ l = . ( ) ( Finite Volume ChPT ) . (60)5 Ansatz c / dofNLO 0.72(46)NLO-fv 1.07(47)Analytic 0.60(44)TABLE XXVIII: Fit ansatze, mass ranges and uncorrelated c / dof obtained in our analyses. The fits wereperformed for pion masses less than 420 MeV.Parameter Value Parameter Value C m p -0.001(1) GeV C m K C m p C m K C m p C f K C f p C f K C f p C f K C f p C f K C f p C m W C m K C m W In table XXIX we present the parameters of the fit with the analytic ansatz over the same massrange as for the fits using SU(2) chiral perturbation theory, as explained in the previous paragraph.We find that analytic fits including a larger range of pseudoscalar masses give an acceptable un-correlated c / dof but then the lightest data points were consistently missed by the fit by aboutone standard deviation. The utility of such extended fits for extrapolating to the physical pointwas therefore compromised and we therefore decided to restrict the range of masses used in theanalytic fits.The global fit to many ensembles of partially quenched data is naturally a high dimensional spaceand so the exposition of the fits is best performed by looking at portions of the data in turn. Inorder to illustrate the quality of the fits, in the following subsections we display the fit and data foreach physical quantity in turn. In total we have analysed five ensembles at two lattice spacings,and each ensemble has measurements at many partially quenched valence-quark masses. As it is6only feasible to present a subset of possible plots, in the following we display the dependence ofeach quantity on the valence quark masses at the lightest sea-quark mass ( m l = .
005 for the 24 ensembles and m l = .
004 on the 32 ensembles). The exception of course, is the mass of theOmega baryon m hhh which does not depend on the light valence-quark masses. We also displaythe unitary subset of data on both lattice spacings along with the mass dependence we infer fromour fits in the unitary continuum limit.Before discussing the chiral and continuum behaviour of hadronic masses and decay constants indetail, we present in table XXX our results for the unrenormalised physical quark masses and thelattice spacings obtained from the three fits. In this table the quark masses are given in latticeunits. The non-perturbative renormalization of the masses will be discussed in Sec. VI where thevalues of the renormalized quark masses in the MS scheme will be presented. NLO NLO fv Analytic˜ m l ( ) m s ( ) a − ( ) . ( ) GeV 2 . ( ) GeV 2 . ( ) GeV˜ m l ( ) m s ( ) a − ( ) . ( ) GeV 1 . ( ) GeV 1 . ( ) GeVTABLE XXX: Unrenormalised physical quark masses in lattice units and the values of the inverse latticespacing a − for the 32 and 24 ensembles.
1. Chiral and continuum behaviour of the W -baryon The W mass is fitted using Eq. (45) (or equivalently (55) ). The fit form for the W baryon doesnot change between the different ans¨atze and only very small differences arise from the differentestimates of physical quark masses and hence of the lattice spacings. For illustration, Figure 27shows the extrapolation of the W mass using the analytic ansatz.7 m l (GeV) m hhh ( G e V ) data24 data FIG. 27: The fit to the light-quark mass behaviour of the W -baryon in the continuum limit obtained usingthe analytic ansatz. The corresponding plots using the infinite and finite-volume SU(2) ChPT ansatz arealmost indistinguishable, differing only slightly in the estimates of the physical quark masses and the latticespacings.
2. Chiral and continuum behaviour of the pion mass
We display the fits of the partially quenched pion masses using infinite volume NLO SU(2) par-tially quenched ChPT (i.e. to the partially quenched generalization of Eq. (38) given in Eq. (B.32)of ref. [1]) in figure 28 for the lightest 24 and 32 ensembles. As discussed in section V C, wedivide by the average valence-quark mass with the intention of enhancing the visibility of chirallogarithms. Figure 29 displays the corresponding fit of the same data but including finite-volumecorrections.It is apparent that the infinite volume and finite volume NLO fits diverge rapidly from our data atlarger masses, and this indeed is the reason why we were compelled to introduce the upper cut-offof 420 MeV for this analysis [1].We now consider the chiral extrapolation of the pion mass using the analytic form of Eq. (49) whichis shown in Fig. 30. Comparing Figs. 28 and 29 with Fig. 30 suggests that data at substantiallylarger masses can be described by the analytic expansion, without any curvature terms in theansatz. The division by the average valence quark mass in the plots, coupled to allowing thetangent not to pass through the origin (i.e. that the extrapolated m p at m x = m y = C m p = m ~ x ( a m xy ) / m ~ a vg m y = 0.04m y = 0.03m y = 0.02m y = 0.01m y = 0.005m y = 0.001 m ~ x ( a m xy ) / m ~ a vg m y = 0.03 m y = 0.025m y = 0.008 m y = 0.006 m y = 0.004 m y = 0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 m ~ x ( a m xy ) / m ~ a vg m y = 0.01m y = 0.005m y = 0.001 m ~ x ( a m xy ) / m ~ a vg m y = 0.008m y = 0.006m y = 0.004m y = 0.002 FIG. 28: Global fits obtained using infinite volume NLO SU(2) chiral perturbation theory for the pion mass.The top-left panel includes the partially quenched data from the m l = .
005 ensemble on the 24 lattice andthe data points in the top-right panel are from the m l = .
004 ensemble from the 32 lattice. In each casethe curves correspond to the appropriate value of the lattice spacing. The points marked by the circles wereincluded in the fit, whereas those marked by the diamonds were not. In the bottom two panels we zoom intothe low-mass region, illustrating the fits to the points which were included (24 points on the left and 32 points on the right). (For fixed ˜ m x , m y decreases as ( am xy ) / ˜ m avg increases.) figure 30 in the unitary chiral limit. In fact we find that C m p is numerically small and consistentwith zero, C m p = − . ( ) GeV . We stress again that while Goldstone’s theorem implies thevanishing of the pion mass in the SU(2) chiral limit, this does not necessarily imply that C m p = C m p = C m p is consistent with zero. This is illustrated by the flat behaviour (within9 m ~ x ( a m xy ) / m ~ a vg m y = 0.01m y = 0.005m y = 0.001 m ~ x ( a m xy ) / m ~ a vg m y = 0.008m y = 0.006m y = 0.004m y = 0.002 FIG. 29: Global fits for the pion mass obtained using NLO SU(2) chiral perturbation theory with finite-volume corrections. In this case we only include the points which were included in the fit ( m l = . points on the left and m l = . points on the right) since the finite-volume corrections at largermasses are small. (For fixed ˜ m x , m y decreases as ( am xy ) / ˜ m avg increases.) the statistical precision) for the chiral behaviour of the unitary points for m p / m l in the continuumlimit shown in the right panel in Fig. 31. Allowing for a non-zero value of C m p does howeverlead to an amplified error for m p / m l at the physical point. The left panel of Fig. 31 shows thecorresponding plots for the infinite and finite-volume ChPT fits.Goldstone’s theorem equally applies at vanishing valence-quark mass ( m x = m y =
0) but with anon-zero sea-quark mass ( m l > C m p was consistent with zero, in the partially quenched direction we find that thecorresponding constant C m p + C m p m l is non-zero, specifically C m p = . ( ) GeV. This value for C m p is much larger than might be created by propagating the mass dependence in m ′ res ( m ) throughthe term involving C m p ; the greatest mass dependence in m ′ res occurs on our 24 ensembles in thepartially quenched direction, but can at most generate a 1% correction to ˜ m and produces a termmuch smaller than the measured C m p . Further, the residual chiral symmetry breaking is four timessmaller for the 32 ensemble which is also included in the global fit. Our results from this globalanalytic fit therefore require a curvature, most likely from partially-quenched chiral logarithmswhich are known to be larger than in the unitary direction, in order for Goldstone’s theorem to besatisfied.It is also worth emphasizing that the discovery of chiral logarithms in lattice data from plots suchas those in Figs. 28 to 30 is to a certain extent artificial. Inconsistency with LO chiral perturbation0 m ~ x ( a m xy ) / m ~ a vg m y = 0.04m y = 0.03m y = 0.02m y = 0.01m y = 0.005m y = 0.001 m ~ x ( a m xy ) / m ~ a vg m y = 0.03 m y = 0.025m y = 0.008 m y = 0.006 m y = 0.004 m y = 0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 m ~ x ( a m xy ) / m ~ a vg m y = 0.01m y = 0.005m y = 0.001 m ~ x ( a m xy ) / m ~ a vg m y = 0.008m y = 0.006m y = 0.004m y = 0.002 FIG. 30: Global fit curves obtained using the analytic fit ansatz (49) overlaying the simulated pion masseson the m l = . ensemble (top-left) and the m l = . ensemble (top-right). Points marked bycircles were included in the fit, those marked by diamonds were not. The simple linear expansion replicatesthe entire range of lattice data reasonably well with the description being rather better than NLO chiralperturbation theory at our larger masses. In the bottom two panels we zoom into the low-mass region,illustrating the fits to the points which were included (24 points on the left and 32 points on the right).(For fixed ˜ m x , m y decreases as ( am xy ) / ˜ m avg increases.) theory is certainly indicated. Our linear fits suggest that the transformations made in displayingthe data render even conclusions of genuine curvature, let alone unambiguous demonstration oflogarithmic mass dependence, to be somewhat optimistic. In order to prove logarithmic behaviour,one should really change quark masses substantially on a logarithmic scale; our present lattice datasupports only the weaker claim of consistency with logarithmic behaviour in the partially quencheddirection.1 m l (GeV) m ll / m l ( G e V ) NLO SU(2) ChPT+FVNLO SU(2) ChPT32 data (ChPT+FV)24 data (ChPT+FV)32 data (ChPT)24 data (ChPT) 0 0.005 0.01 0.015 0.02 m l (GeV) m ll / m l ( G e V ) NLO SU(2) analyticNLO SU(2) ChPT32 data (analytic)24 data (analytic)32 data (ChPT)24 data (ChPT) FIG. 31: Left panel: Pion mass fit for the SU(2) NLO fit form in the continuum limit, both with and withoutfinite volume logarithms. We adjust the data points to the continuum limit using the a dependence in ourfit form and overlay these. Right panel: Chiral extrapolation of the pion mass using the analytic (52) andinfinite-volume NLO ChPT ans¨atze.
3. Chiral and continuum behaviour of the pion decay constant
We now turn to the chiral behaviour of f p and the extrapolation to the physical point. The leadingterm in all the fits contains an a correction and we display the fits performed at non-zero latticespacing combined with the unmodified lattice data and also our continuum predictions combinedwith the lattice data extrapolated to the continuum limit using the results of the fits.We display our fits obtained using infinite volume NLO SU(2) partially-quenched ChPT in Fig-ure 32. The corresponding fits including finite-volume corrections are shown in Figure 33. FinallyFigure 34 displays the fits obtained using our analytic ansatz. Having performed the fits, we adjustour unitary data to the continuum limit using the fitting functions with the determined parametersand display the adjusted data in Fig. 35 together with the finite and infinite-volume NLO SU(2)ChPT fits (left panel) and the analytic fit (right panel). The effect of the adjustment to the con-tinuum limit is illustrated in Figure 36 where the fits are superimposed on the unadjusted unitarydata. It can be seen from Figs. 35 and 36 that the adjustment to the continuum limit for the piondecay constant is very small.The predictions for f p extrapolated to the physical quark masses for each of the fits is given intable XXXI. We anticipate the discussion of the global fits for f K which are presented in Sec V E 6and mention that the predictions for f K extrapolated to the physical quark masses are given intable XXXII, and the predictions for f K / f p extrapolated to the physical quark masses are given in2 m ~ x a f xy m y = 0.04m y = 0.03m y = 0.02m y = 0.01m y = 0.005m y = 0.001 m ~ x a f xy m y = 0.03m y = 0.025m y = 0.008m y = 0.006m y = 0.004m y = 0.0020 0.002 0.004 0.006 0.008 0.01 0.012 0.014 m ~ x a f xy m y = 0.01m y = 0.005m y = 0.001 m ~ x a f xy m y = 0.008m y = 0.006m y = 0.004m y = 0.002 FIG. 32: Global fits to the lattice data for the pion decay constant obtained using infinite-volume NLOSU(2) chiral perturbation theory. The top-left and top-right panels correspond to the 24 , m l = .
005 and32 , m l = .
004 ensembles respectively. Points marked by circles are included in the fits, while those withheavier masses marked by diamonds are not. In the bottom two panels we zoom into the low-mass region,illustrating the fits to the points which were included (24 points on the left and 32 points on the right).(For fixed ˜ m x , m y increases as a f xy increases.) table XXXIII.We find that the NLO SU(2) fits underestimate the physical value at our simulated lattice spacings,and that this discrepancy is amplified a little by the extrapolation to the continuum limit. At eachof our two lattice spacings, the analytic ansatz extrapolates close to the physical value of f p , but,with our ansatz for the form of the a effects, the result becomes statistically inconsistent in thecontinuum limit.From the above discussion we see that using NLO ChPT to perform the chiral extrapolation for f p results in a value which is significantly smaller than the physical one. We recall that only data3 m ~ x a f xy m y = 0.01m y = 0.005m y = 0.001 m ~ x a f xy m y = 0.008m y = 0.006m y = 0.004m y = 0.002 FIG. 33: Global fits to the lattice data for the pion decay constant obtained using NLO SU(2) chiralperturbation theory with finite-volume corrections. In this case we only include the points which wereincluded in the fit ( m l = . points on the left and m l = . points on the right) since thefinite-volume corrections at larger masses are small. (For fixed ˜ m x , m y increases as a f xy increases.)NLO NLO fv Analytic f p f p f continuum p f p in GeV for each global fit ansatz at each simulated lattice spacing and inthe continuum limit. limited to m p <
420 MeV was used in the analysis and note that the fits were performed using thechiral expansion with f , the decay constant in the SU(2) chiral limit, included in the expansionparameter c l / ( p f ) . The downward curvature at low masses seen in Figure 35 can, of course, bereduced by replacing the mass-independent f by an artificial larger parameter such as the physical NLO NLO fv Analytic f K f K f continuum K f K in GeV for each global fit ansatz at each simulated lattice spacing andin the continuum limit. m ~ x a f xy m y = 0.04m y = 0.03m y = 0.02m y = 0.01m y = 0.005m y = 0.001 m ~ x a f xy m y = 0.03m y = 0.025m y = 0.008m y = 0.006m y = 0.004m y = 0.0020 0.002 0.004 0.006 0.008 0.01 0.012 0.014 m ~ x a f xy m y = 0.01m y = 0.005m y = 0.001 m ~ x a f xy m y = 0.008m y = 0.006m y = 0.004m y = 0.002 FIG. 34: Global fits to the lattice data for the pion decay constant obtained using the analytic ansatz inEq. (52). The top-left and top-right panels correspond to the 24 , m l = .
005 and 32 , m l = .
004 ensemblesrespectively. Points marked by circles are included in the fits, while those with heavier masses marked bydiamonds are not. In the bottom two panels we zoom into the low-mass region, illustrating the fits to thepoints which were included (24 points on the left and 32 points on the right). (For fixed ˜ m x , m y increasesas a f xy increases.) f p or f ll ( e m l ) measured at each quark mass used in the simulation. The curvature can also bepartially absorbed by using a subset of terms that arise at NNLO. We have experimented withNNLO fits [46] but find that the low-energy constants are insufficiently constrained by our data tobe of practical use. Thus the resulting predictions for the physical value of f p depend strongly onthe model assumptions used at NNLO.The observed O ( ) deviation found using NLO chiral perturbation theory is broadly consistentwith the size of NNLO terms one might expect to be present at masses in the region of our data.Our data for f p vary from about 20% to 40% above the value of f obtained from our extrapolations5 m l (GeV) f ll ( G e V ) NLO SU(2) ChPT + FVNLO SU(2) ChPT32 data (ChPT+FV)24 data (ChPT+FV)32 data (ChPT)24 data (ChPT) PDG value m l (GeV) f ll ( G e V ) LO analyticNLO SU(2) ChPT32 data (analytic)24 data (analytic)32 data (ChPT)24 data (ChPT) PDG value
FIG. 35: Unitary data for f p adjusted to the continuum limit using each of the fit ans¨atze. The left panelcompares the infinite volume and finite volume forms of the NLO SU(2) fit, while the right panel com-pares the analytic fit to the infinite volume NLO SU(2) fit. The horizontal solid line indicates the value f p − =130.4 MeV (the authors of ref. [45] quote f p − = ( . ± . ± . ) MeV).NLO NLO fv Analytic ( f K / f p ) ( f K / f p ) ( f K / f p ) continuum f K / f p for each global fit ansatz at each simulated lattice spacing and inthe continuum limit. and the square of these terms can be taken as being indicative of the expected NNLO terms. Wemight therefore expect them to be around 5-15% within our simulated mass range.The discrepancy of the prediction for the physical value of f p from the analytic fits is smaller thanthat found with NLO ChPT, but is nevertheless visible. The results at each of the two lattice spac-ings are statistically consistent with f p but lead to an underestimate in the continuum limit. Giventhe sign of the chiral logarithms at NLO, one might expect a linear ansatz to over-estimate ratherthan underestimate the prediction for the physical value. It is nevertheless striking that one cannotadmit any significant non-linearity in this extrapolation and retain consistency with the physicalvalue for f p . The simple analytic form used here appears to be a successful phenomenologicalmodel which is simpler and has fewer parameters than approaches based on ChPT with arbitrarilychosen analytic subsets of NNLO and NNNLO terms.6 m l (GeV) f ll ( G e V ) LO analyticNLO SU(2) ChPTLO analytic 32 LO analytic 24 NLO SU(2) ChPT 32 NLO SU(2) ChPT 24 data24 data PDG value
FIG. 36: Chiral extrapolation of the pion decay constant using the analytic (52) and ChPT (42) fit ans¨atze.Here, the lattice results from the 24 and 32 ensembles are shown along with the mass dependence weinfer both at each lattice spacing and in the continuum limit. The consistency of the two ensembles witheach other and with this continuum limit is indicative of the size of lattice artefacts. The horizontal solidline indicates the value f p − = ( . ± . ± . ) MeV [45].
It is of interest to pose the scientific question whether any of the fit ans¨atze could in principal beconsistent with the experimentally measured pion decay constant? To answer this question weupdate the analysis of Ref. [47] and include an artificially created data point for each ensemblethat represents the experimental result in the continuum limit but includes our fitted a correctionat each non-zero lattice spacing. This is displayed in figure 37 and we find that the analyticans¨atze could be consistent with an uncorrelated c / dof = . ( ) , while NLO ChPT would failto simultaneously fit our data and the physical point, with c / dof = ( ) (infinite volume) and c / dof = ( ) (finite volume).Of course, improved statistical errors, simulations at a third lattice spacing and larger physicalvolumes would give us better control of the continuum extrapolation and finite-volume effects.However, our main conclusion is that it is imperative to simulate with masses substantially nearerto the physical point; this will constrain both fit forms to give more consistent predictions. Ul-timately simulations will be performed directly at physical quark masses and will eliminate thiserror completely. We are currently generating new ensembles with a coarser lattice spacing, with asubstantially larger volume and with very much lighter pion masses (for a preliminary discussion7 m l (GeV) f ll ( G e V ) PDG value m l (GeV) f ll ( G e V ) PDG value
FIG. 37: An artificial data point (the left-most data point in each panel) corresponding to the physical valueof f p [45], but including our uncertainties in the lattice spacing, is added to the data for the pion decayconstant from the five ensembles. The left-hand panel corresponds to the NLO SU(2) ChPT fits and theright-hand panel to the analytic ansatz. of these configurations see Ref. [48]) precisely to address this issue.As an estimate of the systematic uncertainties in physical quantities we take the difference be-tween the results obtained using linear and finite-volume NLO ChPT analyses. This allows for thepossible validity of the full NLO non-analyticity in the region of masses between the data and thephysical point but also recognises that part of this extrapolation may be outside the range of valid-ity of NLO ChPT as suggested by the observation that the present data is surprisingly consistentwith linear behaviour. Guided by the results for f p discussed above, we take as our central valuesfor phenomenological predictions the average of the results obtained from our finite-volume NLOChPT fits and our analytic fits.
4. Chiral and continuum behaviour of the mass of the kaon
We display our fits using infinite volume NLO SU(2) partially quenched ChPT in figure 38. Fig-ure 39 displays the corresponding fits of the same data with the finite-volume corrections included,while the analytic fits are displayed in figure 40. The corresponding unitary view of the data in thecontinuum limit is shown in figure 41. All these plots are for results at the physical sea strangequark mass.8 m ~ x ( a m xh ) m ~ x ( a m xh ) FIG. 38: Dependence of the kaon mass on the mass of the light valence quark with fits performed usinginfinite-volume NLO partially-quenched ChPT. The left panel shows the results from the 24 , m l = . , m l = .
004 ensemble. In each case the results are for the physicalstrange-quark mass. m ~ x ( a m xh ) m ~ x ( a m xh ) FIG. 39: Dependence of the kaon mass on the mass of the light valence quark with fits performed usingfinite-volume NLO partially-quenched ChPT. The left panel shows the results from the 24 , m l = . , m l = .
004 ensemble. In each case the results are for thephysical strange-quark mass.
5. Chiral and continuum behaviour of f K We next discuss f K , the decay constant of the kaon. We display our fits using infinite-volumeNLO SU(2) partially quenched ChPT in Figure 42. The following two figures display fits of thesame partially quenched data to ChPT with finite-volume corrections (Figure 43) and to the globalanalytic fit ansatz (Figure 44). The NLO ChPT fit ans¨atze, both with and without finite-volumelogarithms, are displayed for the unitary data adjusted to the continuum limit in figure 45.9 m ~ x ( a m xh ) m ~ x ( a m xh ) FIG. 40: Dependence of the kaon mass on the mass of the light valence quark with fits performed using theanalytic fit ansatz. The left panel shows the results from the 24 , m l = .
005 ensemble and the right panelfrom the 32 , m l = .
004 ensemble. In each case the results are for the physical strange quark mass. m l (GeV) m l h2 ( G e V ) NLO SU(2) ChPT+FVNLO SU(2) ChPT32 data (ChPT+FV)24 data (ChPT+FV)32 data (ChPT)24 data (ChPT) m l (GeV) m l h2 ( G e V ) LO analyticNLO SU(2) ChPT32 data (analytic)24 data (analytic)32 data (ChPT)24 data (ChPT) FIG. 41: Chiral extrapolation of the kaon mass using unitary data points adjusted to the continuum limit bythe fitting ans¨atze. Here we compare results obtained using the infinite-volume NLO ChPT ansatz to thatusing finite volume logarithms (left panel) and to the analytic ansatz (right panel).
The two panels in Figure 46 display the chiral behaviour of the actual unitary data from the twosets of ensembles (left panel) as well as of the data adjusted to the continuum limit (right panel).From these fits our final predictions for f K are given in table XXXII, and the corresponding resultsfor f K f p in table XXXIII.
6. Predictions
We now present our results for f p , f K and their ratio as well as for the physical bare quark masses.As discussed above, our central value for any physical quantity is taken to be the average of the0 m ~ x a f xh m ~ x a f xh FIG. 42: Dependence of the kaon decay constant on the mass of the light valence quark with fits performedusing infinite-volume partially quenched NLO ChPT. The left panel shows the results from the 24 , m l = .
005 ensemble and the right panel from the 32 , m l = .
004 ensemble. In each case the results are for thephysical strange quark mass. results obtained from analyses using the NLO SU(2) ChPT fit with finite volume corrections andthose from the analytic fit. The difference between the analytic and finite-volume NLO SU(2) fitsis taken as a systematic error. This procedure includes a NLO finite-volume correction, estimatedfrom the difference between results obtained using NLO ChPT at infinite and finite volumes, andwhich is much smaller than the total systematic error here.Our predictions for pseudoscalar decay constants therefore contain systematic errors for finitevolume effects, the chiral extrapolation, and residual chiral symmetry breaking, while the discreti-sation error is included indirectly by the fitting procedure: f continuum p = ( )( ) MeV (61) f continuum K = ( )( ) MeV (62) ( f K / f p ) continuum = . ( )( ) , (63)where we display the statistical and systematic errors separately. We note that the known, exper-imental value of f p influenced our choice to take the central value of physical quantities as theaverage of the results from the analytic and finite-volume NLO ChPT ans¨atze. The prediction for f p cannot therefore be considered unbiased, however as our aim is to select the most likely centralvalue for phenomenologically important quantities such as f K / f p and B K our procedure is bothappropriate and contains a prudent systematic error.1 m ~ x a f xh m ~ x a f xh FIG. 43: Dependence of the kaon decay constant on the mass of the light valence quark. The left panelshows the results from the 24 , m l = .
005 ensemble and the right panel from the 32 , m l = .
004 ensemble.In each case the results are for the physical strange quark mass. There are two curves plotted. The orangecurve is the result one infers for the infinite volume, while the red curve is the result we obtain on the finitevolume. As we do not adjust our data for finite volume effects, the red curve should go through our data.The orange curve also goes through our data which is an indication that the finite volume effects in ourdata are substatistical, and the difference between the orange and red curves at lighter masses indicates thatone should expect substantial finite volume effects if one were to simulate at these lighter masses withoutchanging our present volume. Applying the same procedure to obtain predictions for the physical bare quark masses for the b = .
25 32 ensembles, we find:˜ m ud = . ( )( ) MeV and ˜ m s = . ( )( ) MeV , (64)and these will be renormalised in the following section. The corresponding bare masses for the b = .
13 24 ensembles can be obtained by dividing the results in (64) by the values of Z l and Z h in Table XXVI.
7. Chiral and continuum behaviour of r and r Finally in this section we apply the combined chiral/continuum extrapolation procedure to thescales r and r . Assuming a linear dependence for the light sea-quark mass dependence, andincluding a leading order a term as before, the scales are independently fit to the form r i = c r i + c r i , a a + c r i , m l ˜ m l , (65)2 m ~ x a f xh m ~ x a f xh FIG. 44: Dependence of the kaon decay constant on the mass of the light valence quark with fits performedusing the analytic fit ansatz. The left panel shows the results from the 24 , m l = .
005 ensemble and theright panel from the 32 , m l = .
004 ensemble. In each case the results are for the physical strange quarkmass. m l (GeV) f l h ( G e V ) NLO SU(2) ChPT+FVNLO SU(2) ChPT32 data (ChPT+FV)24 data (ChPT+FV)32 data (ChPT)24 data (ChPT) FIG. 45: Chiral extrapolation of the kaon decay constant for unitary data in the continuum limit. Wecompare the NLO ChPT ansatz to the corresponding ansatz with finite-volume logarithms. where i = ,
1. Prior to the fit, the data are linearly interpolated to each of the physical strangequark masses obtained from the global fits and presented in Table XXX, and the fit and the subse-quent extrapolation are performed using the corresponding physical light-quark mass and latticespacings.The parameters and c / d . o . f of the fits are given in Tables XXXIV and XXXV respectively, andplots showing the fits overlaying the data in the continuum limit are shown in figure 47. The fits3 m l (GeV) f l h ( G e V ) LO analyticNLO SU(2) ChPT32 data (analytic)24 data (analytic)32 data (ChPT)24 data (ChPT) 0 0.005 0.01 0.015 0.02 m l (GeV) f l h ( G e V ) LO analyticNLO SU(2) ChPT32 data (analytic)24 data (analytic)32 data (ChPT)24 data (ChPT) FIG. 46: Chiral extrapolation of the kaon decay constant for unitary data in the continuum limit. Wecompare the NLO ChPT ansatz to the analytic ansatz. The left panel displays the data and fits at non-zerolattice spacing, while the right panel displays the predicted results and correspondingly adjusted data pointsfor the continuum limit. (a) r Parameter ChPT ChPT-fv Analytic c r − − − c r , a -0.25(14) GeV -0.25(14) GeV -0.25(14) GeV c r , m l − − − (b) r Parameter ChPT ChPT-fv Analytic c r − − − c r , a -0.15(11) GeV -0.15(11) GeV -0.15(12) GeV c r , m l -1.76(64) GeV − -1.76(64) GeV − -1.76(64) GeV − TABLE XXXIV: Parameters of the chiral/continuum fits to r and r . to r appear to describe the data well by eye, and have a reasonable (uncorrelated) c / d . o . f forthe central value, but with a large deviation across the superjackknife distribution. The fits to r also appear to describe the data reasonably well, although there does seem to be a tension with theheaviest point on the 24 ensembles, which is likely responsible for the larger c / d . o . f. As thereare only five data points it is difficult to reach any stronger conclusions regarding the data: moreensembles and better statistics are needed. For the purpose of quoting a final result, we apply aPDG scale factor of p c / d . o . f to the statistical errors on each of the results. In order to retain4 Quantity ChPT ChPT-fv Analytic r r c / d . o . f of the chiral/continuum fits to r and r .Quantity ChPT ChPT-fv Analytic r − − − r − − − r / r r and r and the ratio r / r at physical quark masses determinedfrom a chiral/continuum fit using the lattice spacings and quark masses obtained from the global fits. the correlations between these quantities when the ratio is taken, the scale factor is applied to thedifference of each jackknife sample from the mean.The continuum results for r , r and their ratio at physical quark masses are given in table XXXVI.Using the procedure for combining the results obtained using the different chiral ans¨atze outlinedin Section V E 3 and applying the PDG scale factor as above, gives: r = . ( ) stat ( ) FV ( ) c GeV − = . ( ) stat ( ) FV ( ) c fm , r = . ( ) stat ( ) FV ( ) c GeV − = . ( ) stat ( ) FV ( ) c fm , and r / r = . ( ) stat ( ) FV ( ) c , (66)where the finite volume error arising from the different determinations of the lattice spacings andquark masses is smaller than the quoted precision on the ratio. c labels the error due to the chiralextrapolation. For comparison, the MILC collaboration recently obtained r = . ( )( + − ) fm( ≃ . ( )( + − ) GeV − ) [49], and also r = . ( )( ) fm ( ≃ . ( )( ) GeV − ) and r = . ( )( ) fm ( ≃ . ( )( ) GeV − ) from an earlier study [50]. At this time we do not havean explanation of the discrepancy between our results in (66) and those of the MILC collabora-tion beyond noting the very different approaches to setting the scale and performing the chiralextrapolation.5 m l (GeV) r ( G e V ) - data24 data 0 0.005 0.01 0.015 0.02 m l (GeV) r ( G e V ) - data24 data FIG. 47: The scales r (left) and r (right) corrected to the continuum limit, overlaid by the chiral/continuumfit. The extrapolated point at the physical light quark mass is shown as the grey cross. Here the latticespacings and physical light quark mass were obtained from the global fits using the analytic ansatz. Thefits using the quantities obtained with the ChPT and ChPT-fv global fit ans¨atze are almost indistinguishablefrom those shown in these figures. VI. LIGHT-QUARK MASSES
The quark masses quoted in Eq. (64) are the bare masses for the lattice action which we are usingon the 32 ensembles with b = .
25 corresponding to a lattice spacing a − ≃ .
28 GeV. In orderto be useful in phenomenological applications these results must be translated into renormalizedmasses in some standard continuum scheme. Therefore in Subsection VI A we determine therenormalization constants relating the bare masses in (64) to those renormalized in the MS schemeat a renormalization scale of 2 GeV. In Subsection VI B we then combine these renormalizationconstants with the bare masses in (64) to obtain the renormalized masses, the LO LEC B MS ( ) and the chiral condensate. A. Non-perturbative renormalization for quark masses
The quark-mass renormalization factor which relates the lattice bare quark mass to that in the MSscheme is determined using non-perturbative renormalization (NPR) with the RI/SMOM schemesproposed in Ref. [14] as intermediate schemes. This is an extension of the Rome-SouthamptonNPR program in which the RI/MOM scheme was defined [51]. Quark masses renormalized in theRI/SMOM or RI/MOM schemes are obtained entirely non-perturbatively. Since it is not possible to6simulate in a non-integer number of dimensions, continuum perturbation theory is needed to matchthe results in either the RI/SMOM or the RI/MOM scheme and the target MS scheme. We stresshowever, that we completely avoid the use of lattice perturbation theory which often convergesmore slowly than continuum perturbation theory (PT). Since RI/MOM and any of the schemesproposed in [14] are legitimate renormalization schemes, we exploit the freedom to choose anintermediate scheme to reduce its effect on the final result for the renormalized quark mass in theMS scheme and to have a better understanding of this uncertainty.Our earlier study [13], used to normalize the quark mass on the 24 ensembles, applied theRI/MOM scheme to renormalize the quark masses and suffered from sizable systematic errorswith two dominant sources. One of these is the truncation error in the perturbative continuummatching between the RI/MOM and MS schemes. This was estimated to be 6% for m = a − = .
73 GeV and m = Z m is conveniently calculated using the relation Z m = / Z S = / Z P , (67)where Z m , Z S , Z P are the quark mass, flavor non-singlet scalar and pseudoscalar renormalization7factors respectively. Here we are exploiting the important chiral symmetry properties of DWF. Ourconvention is that the renormalization factors multiply the bare quantities to yield renormalizedones: m R = Z m e m , P aR = Z P P a , S aR = Z S S a , (68)where the left-hand sides are the renormalized mass, pseudoscalar and scalar densities and a is aflavour label. e m in Equation (68) is in physical units. The relations in Eq. (67) are necessary for theWard-Takahashi identities to hold for the renormalized operators. The RI/MOM renormalizationcondition on the amputated scalar vertex P S reads Z S Z q
112 Tr [ P S · I ] = . (69) Z q is the wave function renormalization factor, which can be determined using the trace conditionon the local vector operator, Z V Z q
148 Tr [ P V m · g m ] = . (70)The vertex functions P depend on the incoming and outgoing momenta on the two fermion lines, P ( p in , p out ) . The conventional RI/MOM scheme is defined using the forward vertex with p in = p out = p . The renormalization conditions Eqs. (69), (70) are applied by setting the renormalizationscale m to be the off-shell external momentum, m = p , in the chiral limit.It is in principle possible to determine Z S (= Z P ) using the pseudoscalar vertex function instead ofthe scalar one in Eq. (69). However, with the original RI/MOM choice for the external momenta,the pseudoscalar vertex couples to the zero-momentum pion, and the Green function diverges as1 / m q as the quark mass m q → p [54]. Therefore the pseudoscalar vertex cannot beused without some manipulation of the divergence (see e.g. [55]) and has not been considered inour previous publication [13]. This is in contrast with the RI/SMOM schemes described belowwhich do not have such a pole as m q →
0. Similarly, the axial-vector vertex can be used todetermine Z q because Z V = Z A . However, Z q obtained using the vector and axial-vector vertices atlarge but finite p will differ because of the coupling of the axial current to the Goldstone boson[51]. These differences are known to be of O ( / p ) at high momentum from the operator productexpansion [51, 54] or from Weinberg’s theorem of power counting for a Feynman diagram [13].In Ref. [13], the average of the vector and the axial-vector vertex was used to determine Z q andthe difference was included in the systematic error, though the corresponding 1% error is sub-dominant.8The caveats mentioned in the two preceding paragraphs are both connected to the RI/MOMscheme and its channel with an “exceptional momentum”; specifically, the momentum transfer q ≡ p in − p out =
0. This is the reason for the large NPE error. It was demonstrated that the useof non-exceptional momenta p in − p out = p in = p out = q and hence introduce the name “Symmetric Mom”(SMOM) schemes. The two schemes RI/SMOM and RI / SMOM g m are defined with this kine-matical choice but differ in the G -projection operators which are used to define the wave functionrenormalization. For the vector (axial-vector) vertex function the projector q / q m / q ( g q / q m / q ) isused in the RI/SMOM scheme and g m ( g g m ) as in Eq. (70) is used for RI / SMOM g m . The standard I ( g ) spinor projector is used for the scalar (pseudoscalar) vertex in both new schemes.The conversion factors from the RI/SMOM and RI / SMOM g m schemes to MS have been calculatedat one-loop order in Ref. [14] and recently to two-loop order [15, 16]: C m ( RI / SMOM → MS , m ) = − (cid:18) a s ( m ) p (cid:19) . − (cid:18) a s ( m ) p (cid:19) ( . + . n f ) · · · , (71) C m ( RI / SMOM g m → MS , m ) = − (cid:18) a s ( m ) p (cid:19) . − (cid:18) a s ( m ) p (cid:19) ( . + . n f ) · · · , (72)where the coefficients have been rounded to the third decimal place. Evaluating these factors at m = C m ( RI / SMOM → MS , m = , n f = ) = − . − . · · · , (73) C m ( RI / SMOM g m → MS , m = , n f = ) = − . − . · · · . (74)In the RI/MOM and RI ′ /MOM schemes the conversion factors are known to three-loop order [56,57]: C m ( RI / MOM → MS , m = , n f = ) = − . − . − . + · · · , (75) C m ( RI ′ / MOM → MS , m = , n f = ) = − . − . − . + · · · . (76)We note that, at least up to two-loop order, the convergence of the series relating the new SMOMschemes to MS is considerably better than for the RI/MOM scheme. As already mentioned, thetruncation error of the RI/MOM scheme was estimated from the size of the highest order termavailable (3 loop). Having in addition two intermediate SMOM schemes, we can expect to have amore reliable estimate of the truncation error.9We now turn to the numerical evaluation of the renormalization factors. At each value of b , weuse data obtained at the three light-quark masses: m l = . .
006 and 0 .
008 for the finer 32 lattice and m l = . .
01 and 0 .
02 for the coarser 24 lattice. 20 configurations were analyzedfor each point. The ratio of quark wavefunction and local axial current renormalization factors iscalculated from the average of vector and axial-vector vertex functions, Z q Z V = ( L V + L A ) , (77)with projected and traced vertex functions: L RI/SMOM V =
112 ˆ q Tr [ P V m · ˆ q / ˆ q m ] and L RI/SMOM A =
112 ˆ q Tr [ P A m · g ˆ q / ˆ q m ] , (78)for the RI/SMOM scheme. Here q m in the continuum RI/SMOM scheme [14] has been replacedwith the ˆ q m = sin ( q m ) , as the derivative for the divergence of the current in the continuum theoryis naturally replaced by the symmetric difference on the lattice. A remarkable feature of theRI/SMOM scheme is that in the chiral limit L V = L A holds non-perturbatively, in contrast to L V = L A for RI/MOM scheme due to spontaneous symmetry breaking (SSB). In principle there couldstill be a small difference for the lattice RI/SMOM scheme with non-zero m res , which, however, isnegligible in the momentum range we use [52]. Using the continuum Ward-Takahashi identities,one can also show the equivalence of Z q in the RI/SMOM and RI ′ /MOM schemes [14].The RI / SMOM g m scheme is defined using the conventional projectors, L RI/SMOM g m V =
148 Tr [ P V m · g m ] and L RI/SMOM g m A =
148 Tr [ P A m · g g m ] . (79)Although these projectors are superficially the same as those used in the RI/MOM scheme, itshould be remembered that the kinematics is different in the two cases with no exceptional chan-nels in the Green functions used to define the RI / SMOM g m scheme.The product of mass and wavefunction renormalization factors is calculated from the average ofscalar and pseudoscalar vertex functions, Z m Z q = ( L S + L P ) , (80)with L S =
112 Tr [ P S · ] and L P =
112 Tr [ P P · g ] , (81)again defined with the SMOM kinematics for the vertex functions. While L S = L P holds to allorders in perturbation theory with naive dimensional regularization, by using Weinberg’s power-counting scheme we see that they can in general differ by terms of O ( / p ) [13]. The difference0 [GeV ]0.0010.01 L P - L S FIG. 48: L P − L S as a function of p [GeV ] for fine (32 ) and coarse (24 ) lattices. A straight line with1 / p slope but arbitrary normalization is drawn to guide the eye. L P − L S after the chiral extrapolation is plotted in Fig. 48 as a function of p (in physical units)for both the 24 and 32 lattices. The figure confirms the expected approximate 1 / p scaling. Theunwanted non-perturbative effect from SSB is small and the introduction of non-exceptional mo-menta has had the expected effect. This is in contrast to the RI/MOM scheme with the exceptionalchannel, where the same difference behaves as 1 / ( mp ) , and thus diverges in the chiral limit atfinite p .The mass renormalization factor Z s m , with s = RI / SMOM or RI / SMOM g m , is given by combiningEqs. (77) and (80), Z s m = Z V L S + L P L s V + L s A . (82)In calculating the ratio of vertex functions in Eq. (82) we take the average of S and P or V and A for each light-quark mass and then fit with a quadratic ( c + c ′ ( m l + m res ) ) or linear c + c ′′ ( m l + m res ) formula to obtain the value c in the chiral limit for the numerator and denominator. For illustration,the extrapolation for the numerator using the quadratic formula is shown in Fig. 49, where theobserved mass dependence is seen to be very small. Because of the very mild mass dependence,to the precision with which we quote our results and errors, the quadratic and linear extrapolationformulae lead to exactly the same quark-mass renormalization factor and error. Finally taking theratio and combining with Z V gives the mass renormalization factor in the RI/SMOM schemes. Therenormalization factor in the MS scheme at a scale m = s to MS at m = p in = p out = q using Eqs. (71) and (72) and then running to 2 GeV using1 f +m res p /( p /16) = FIG. 49: Chiral extrapolation of ( L P + L S ) / ) lattice for each p point. Z m fi SMOM g m (p) fi SMOM(p) fi SMOM g m fi MS(2GeV) fi SMOM fi MS(2GeV) p = G e V FIG. 50: Z SMOM g m m ( m ) and Z SMOM m ( m ) as functions of m = p , and Z MS m ( ) from the SMOM orSMOM g m schemes as function of matching scale squared p for the fine lattice. The interpolation points areshown with the error bar at p = ( ) . the three-loop anomalous dimension in the MS scheme. We use the four-loop QCD beta functions[58] to calculate a ( ) s ( m ) for running and matching as shown in Appendix A of Ref. [13]. Therelevant parameters taken from the 2008 Particle Data Group [45] are a ( ) s ( m Z ) = . , m Z = . , m b = .
20 GeV and m c = .
27 GeV , (83)where the quark masses are in the MS scheme at the scale of the mass itself, e.g. m b = m MS b ( m b ) .In Fig. 50 we plot Z SMOM g m m ( m ) and Z SMOM m ( m ) in the SU(2) chiral limit as functions of m = p ensembles. In addition we also plot Z MS m ( ) as functions of the matching scale p obtained with SMOM and SMOM g m as the intermediate schemes. In an ideal situation, i.e. one inwhich the errors due to NPE contamination, truncation of perturbation theory and lattice artifactsare all small, the results obtained using the two intermediate schemes would give the same resultsfor Z MS m ( ) , and the results would be independent of ( pa ) . Since we have observed that theNPE error is small, the difference between the two sets of results is mostly due to the truncation ofperturbation theory and lattice discretization errors. The observed decrease in this difference as p increases is consistent with the expected behaviour of the truncation error. Conversely, since thetruncation error increases as p decreases, taking the limit ( pa ) →
0, which is a typical treatmentto eliminate the discretization error, is not an appropriate procedure. We therefore choose insteadto evaluate Z m by taking an intermediate reference point p = ( ) , for both the 24 and 32 lattices. In this way, as we take the continuum limit of the renormalized quark mass, the leading ( pa ) discretization error associated with the non-perturbative renormalization will be removed.There is a subtlety due to lattice artefacts which are not O ( ) invariant and which are responsiblefor the non-smooth ( pa ) dependence in the figure. A term like a (cid:229) m ( p m ) / p , whose presencehas been demonstrated in the conventional RI/MOM scheme for Wilson quarks [59], could existalso in the SMOM schemes. Such a term would manifest itself as scattered data around a smoothcurve in p , and the size of the scatter is expected to be comparable to the leading ( pa ) error asboth are of the same order in a . This appears to be compatible to what is shown in the figure.Of course, it would be very helpful to know these terms, but in the absence of this knowledge weinclude this scatter in the systematic error by inflating the error by a factor p c / dof. The resultsare Z MS ( ) m ( m = , n f = SMOM g m ) = . ( ) , (84) Z MS ( ) m ( m = , n f = SMOM ) = . ( ) . (85)The final arguments on the left-hand sides denote the choice of intermediate scheme. The error onthe right-hand sides is the combination of the statistical fluctuations and the scatter of the pointsaround the linear fit. The central values and errors are shown in the figure at the reference point, p = ( ) .The 24 coarser lattice has been analyzed similarly for the m l = . .
01 and 0 .
02 ensemblesand the results are shown in Fig. 51. The mass renormalization factors on the 24 lattice for the3 Z m p = G e V FIG. 51: Same figure as Fig. 50, but for the coarse 24 lattice. two intermediate SMOM schemes are: Z MS ( ) m ( m = , n f = SMOM g m ) = . ( ) , (86) Z MS ( ) m ( m = , n f = SMOM ) = . ( ) . (87)In Eq. (64) we have presented the bare quark masses for the fine 32 lattice and in Table XXVIwe give the ratios of equivalent bare masses on the 24 and 32 lattices. Because of the different O ( a ) artefacts for the light and heavy quark masses, there are two such ratios Z l for the ud quarksand Z h for the s quark. These ratios Z l and Z h are also the scheme-independent ratios of therenormalization constants on the course and fine lattices. We now use these ratios to estimatethe difference of the MS renormalized masses with the SMOM and SMOM g m schemes in thecontinuum limit. The continuum extrapolation of Z ( ) m and Z ( ) m / Z l or Z ( ) m / Z h will remove the ( pa ) error in the non-perturbative renormalization. Thus, if a difference is found, it can largely beattributed to the truncation error of the perturbative matching. Performing such an extrapolationwe find Z MS ( ) cml ( m = , n f = SMOM g m ) = . ( ) , (88) Z MS ( ) cml ( m = , n f = SMOM ) = . ( ) , (89)for the ud quark, and Z MS ( ) cmh ( m = , n f = SMOM g m ) = . ( ) , (90) Z MS ( ) cmh ( m = , n f = SMOM ) = . ( ) (91)4 ensemble fine (32 ) coarse (24 ) coarse (16 )[13]intermediate scheme RI / SMOM RI / SMOM RI / MOMPT truncation error 2.1% 2.1% 6% m s = ( L P − L S ) / ¥ ) ( L A − L V ) / Z MS m ( ) with intermediate RI / SMOM schemes (thiswork) and RI / MOM scheme [13]. for the s quark. Note that because these factors multiply ˜ m ud ( ) / a ( ) or ˜ m s ( ) / a ( ) presented in Eq. (64) to give the MS mass in the continuum limit, they are made to absorb the O ( a ( )) discretization error in these bare quark masses on the fine lattice. Because of this,as well as the fact that the Z m ’s are free from O ( a ) errors originating from the SMOM non-perturbative renormalization, we have put additional suffix “ c ” as “continuum” to distinguish themfrom Z MS ( ) m . The existence of a mass dependent contribution to the O ( a ) artefacts gives rise tothe different Z m for the light and heavy-quark masses. From the two different estimates of the MSrenormalization factors with the SMOM and SMOM g m intermediate non-perturbative schemes,we choose to take SMOM g m for our central value. The reason is that the scatter about the linearbehaviour observed for the SMOM scheme in Figs. 50 and 51 is much larger. Although the effectof the scatter has been taken into account in the error, we consider the continuum extrapolationfrom the SMOM scheme to be less reliable. The difference in the central values of Z MS ( ) cml inEqs. (88) and (89) is about 1%, and this is also the case for the difference between the centralvalues of Z MS ( ) cmh in Eqs. (90) and (91). These differences of about 1% give an indication of thepossible size of the truncation error of the perturbative two-loop matching to MS (it should benoted however, that the errors in the renormalization factors in the SMOM scheme are even alittle larger). Another estimate of the truncation error of the matching is obtained by evaluatingthe size of the two-loop term in Eq. (74), resulting in 2 .
1% for the SMOM g m scheme. In orderto be conservative, we shall take the latter as our estimate. Other systematic errors arise from thefact that the simulated strange mass is non-zero and from the small difference in the scalar andpseudoscalar vertices due to the residual spontaneous symmetry breaking effects. The first error5is estimated from the response of scalar and pseudoscalar vertex functions to the variation of thelight-quark mass [13]. From the flat behaviour of L P + L S on the light-quark mass in Fig. 49 it canbe seen that this uncertainty is small. The error estimates are compiled in Table XXXVII. In thetable, the corresponding errors from the RI / MOM analysis [13] are shown for comparison. Allerrors have become significantly smaller for the new SMOM schemes. Now our final values forthe MS renormalization factor read Z MS ( ) cml ( m = , n f = ) = . ( )( ) , (92) Z MS ( ) cmh ( m = , n f = ) = . ( )( ) , (93)where the first error is the statistical uncertainty inflated to take into account the scatter about thelinear behaviour due to O ( ) non-invariant effects (as explained above) and the second is due tothe remaining systematic effects and is dominated by the 2 .
1% truncation error of the perturba-tive matching. Here we have not taken into account the statistical fluctuation of Z V , which will beproperly included in the calculation of the renormalized quark masses described in the next subsec-tion. The corresponding renormalization factor for the light-quark mass on the coarse 24 latticeis Z MS ( ) cml ( m = , n f = ) = Z l · Z MS ( ) cml ( m = , n f = ) = . ( )( ) . This valueis consistent with our earlier estimate of the same quantity using RI/MOM as the intermediatescheme, 1 . ( ) [13], but now with a considerably reduced error. B. Renormalized quark masses
After the detailed discussion of the quark-mass renormalization, it is now straightforward to com-bine the renormalization constants in Eqs. (92) and (93) with the physical bare quark masses onthe 32 lattice in Eq. (64) to obtain the light and strange quark masses renormalized in MS scheme: m MS ud ( ) = Z MS ( ) cml ( m = , n f = ) · ˜ m ud ( ) · a − ( )= . ( ) stat ( ) sys ( ) ren MeV , (94) m MS s ( ) = Z MS ( ) cmh ( m = , n f = ) · ˜ m s ( ) · a − ( )= . ( . ) stat ( . ) sys ( . ) ren MeV , (95)where the three errors on the right-hand side correspond to the statistical uncertainty, the system-atic uncertainty due to the chiral extrapolation and finite volume, and the error in the renormaliza-tion factor. We recall that for the error due to the chiral extrapolation we conservatively take the6full difference of the results obtained using the finite-volume NLO SU(2) and analytic fits and forthe central value we take the average of these results. We estimate the finite-volume effects fromthe difference of the results obtained using finite volume and infinite-volume NLO ChPT fits andcombine these errors in quadrature. The finite-volume errors prove to be small. The error in therenormalization factor includes those in Eqs. (92) and (93).The ratio of the s and ud quark masses is m s m ud = . ( . ) stat ( . ) sys . (96)We end this section by presenting our results for the leading-order LEC B and the chiral conden-sate. Using the finite-volume NLO ChPT fits we find B MS ( ) = Z MS ( ) − ml ( m = , n f = ) · B ( ) · a − ( ) = . ( ) stat ( ) sys ( ) ren GeV . (97)Combining this result with the pion decay constant in the chiral limit, also obtained using thefinite-volume NLO ChPT fits the chiral condensate is found to be [ S MS ( )] / = [ f B ( ) / ] / = ( ) stat ( ) sys ( ) ren MeV . (98)In Eqs. (97) and (98) the second error is only due to finite volume corrections estimated from thedifference of finite and infinite volume NLO ChPT fits. VII. TOPOLOGICAL SUSCEPTIBILITY
The topological charge Q , defined on a single Euclidean space-time configuration, and its sus-ceptibility, c Q , are interesting quantities to calculate. While Q depends only indirectly on thequark masses, leading order SU(2) ChPT [60, 61] predicts a strong dependence of c Q on the lightsea quark mass with c Q vanishing linearly as m l →
0, suggesting that c Q may show importantdynamical quark mass effects.In the continuum Q and c Q are defined by Q = g p Z d x G mn ( x ) ˜ G mn ( x ) and c Q = h Q i / V , (99)where V is the four-volume of the lattice, G mn ( x ) is the gluon field strength tensor and ˜ G mn ( x ) , itsdual. In the continuum, Q is integer valued and related to exact chiral zero modes of the masslessDirac operator by the Atiyah-Singer index theorem [62]. For sufficiently smooth gauge fields it7is possible to find a lattice expression which will always evaluate to an integer [63], as in thecontinuum limit. However, in the calculation reported here the necessary smoothness conditionis not obeyed and we instead replace the right-hand side of Eq. (99) by a sum of Wilson loopschosen to approximate the G mn ( x ) ˜ G mn ( x ) product in Eq. (99). Specifically we employ the “five-loop improved” (5Li) definition of the topological charge proposed in Ref. [64] which at tree levelis accurate through order a . However, before evaluating this lattice expression for the topologicalcharge, we smooth the links in the lattice by performing a series of APE smearing steps [65,66]. The smearing parameter was set to 0.45, and 60 smearing sweeps were performed beforemeasuring Q . The results are insensitive to the choice of these parameters.In Fig. 52 the Monte Carlo time history of Q is shown for each ensemble of gauge fields in ourstudy. For each case, the update algorithm RHMC II [1] was used, except for the first 1455configurations for the m l = .
01 ensemble where the RHMC 0 and RHMC I algorithms wereused. In [1] it was shown that RHMC II is more effective in changing the gauge field topology,and therefore produces shorter auto-correlation times. The data for the first half (up to trajectory5000) of both 24 ensembles is repeated from [1]. Figure 52 shows clearly the expected slowingof the rate of change of topological charge when moving towards the continuum [67] and, to alesser degree, when decreasing the quark mass. The integrated auto-correlation times for Q forthe smaller lattice spacing ensembles are shown in Fig. 2. While this figure is consistent with theautocorrelation times reaching a plateau of about 80 time units when integrated over an intervalof about 200 time units, the exploding errors make this conclusion highly uncertain. ScanningFig. 52 by eye, one might argue that the auto-correlations could be 500 time units, or longer. Forexample, note the large fluctuation to negative Q beginning around time unit 4750 for m l = . | Q | .Because of the parity symmetry of our calculation, the average of the pseudo-scalar quantity h Q i vanishes. However, c Q remains non-zero and at leading order in SU(2) chiral perturbation the-ory [60, 61] is given by c Q = S (cid:18) m u + m d (cid:19) − = S m u m d m u + m d , (100)where S = B f / c Q = S (cid:18) m u + m d (cid:19) − × (cid:18) − ( p f ) m p log m p L + K ( m u + m d ) + ( K + K ) m u m d m u + m d (cid:19) , (101) = S m l (cid:18) − ( p f ) m ll log m ll L + ( K + K + K ) m l (cid:19) , (102)where K i = S L i / f are proportional to the Gasser-Leutwyler NLO LEC’s [68], and in the lastline the formula is evaluated for degenerate quarks. In contrast to other quantities considered inthis paper, we do not attempt to characterize or evaluate the corrections to Eqs. (101) or (102)which come from non-zero lattice spacing. That interesting question is left for future work.In Tab. XXXVIII values of h Q i and c Q for each ensemble of configurations are summarized. Totest for the expected auto-correlations, the data were blocked into bins of various sizes rangingfrom 10 to 600 time units. The quoted values of the statistical errors resulted when the block sizeswere taken large enough that the errors no longer changed significantly. The block sizes are givenin Tab. XXXVIII. For all cases the first 1000 time units were discarded for thermalization.The dependence of c Q on the light quark mass is shown in Fig. 54. All of the data points lieabove the LO curve (dashed line), all but the lightest significantly so. The result of the fit ( c / dof ≈ / ≈
3) to the NLO formula Eq. (102) is also shown. Since we have not determined K in Eq. (102) from other means, we treat the linear combination of LEC’s as a single, new, freeparameter in the fit and find ( K + K + K ) = . ( . ) . Except for the lightest data point,there is scant evidence for large O ( a ) errors, though the statistical errors on the heavier twopoints with a − = .
284 are somewhat large. Omitting the former point in the fit leads to a moreacceptable value of c /dof ≈ .
5, suggesting the lightest point may be systematically low due tolong auto-correlations in Q that are not well resolved in our finite Markov chain of configurations.Despite these limitations, the data appear to show a dependence on the light sea quark mass that isconsistent with the dictates of NLO SU(2) ChPT. VIII. CONCLUSIONS
We have presented results from simulations using DWF and the Iwasaki gauge action for latticeQCD at two values of the lattice spacing ( a − = 1.73 (3) GeV and a − = 2.28 (3) GeV) and for uni-9 TABLE XXXVIII: Topological charge and susceptibility. The measurement frequency, “meas. freq.”, and“block size” are given in units of Monte Carlo time. m l meas. freq. block size h Q i h Q i c (GeV )0.005 5 50 0.49 (25) 28.6 (1.4) 0.000290 (14)0.01 5 50 -0.22 (37) 45.2 (2.5) 0.000458 (25)0.004 4 200 0.59 (42) 11.4 (1.1) 0.000148 (14)0.006 4 200 -0.07 (64) 24.8 (4.3) 0.000322 (55)0.008 4 400 0.64 (100) 27.9 (5.6) 0.000363 (72) -20-1001020-20-1001020-20-1001020-20-10010200 1000 2000 3000 4000 5000 6000 7000 8000 9000-20-1001020 FIG. 52: Monte Carlo time histories of the topological charge. The light sea quark mass increases from topto bottom, (0.005 and 0.01, 24 (top two panels), and 0.004-0.008, 32 ). Data for the 24 ensembles up totrajectory 5000 were reported originally in [1] and the results from the new ensembles are plotted in black.Most of the data was generated using the RHMC II algorithm (red and black lines). The RHMC 0 (greenline) and RHMC I (blue line) algorithms were used for trajectories up to 1455 for the m l = .
01 ensemble.The small gap in the top panel represents missing measurements which are irrelevant since observables arealways calculated starting from trajectory 1000. -20 -15 -10 -5 0 5 10 15 20Q050100150200 -20 -10 0 10 20Q050100150200 -20 -10 0 10 20Q050100150200-20 -15 -10 -5 0 5 10 15 20Q050100150200 -20 -10 0 10 20Q050100150200 FIG. 53: Topological charge distributions. Top: 32 , m l = . − . , m l = .
005 and 0.01. tary pion masses in the range 290–420 MeV (225–420 MeV for the partially quenched pions). Theraw data obtained at each of the two values of b was presented in Sections III and IV respectivelyand the chiral behaviour of physical quantities on the 24 and 32 lattices separately was studiedin Appendix A. The main aim of this paper however, was to combine the data obtained at thetwo values of the lattice spacing into global chiral–continuum fits in order to obtain results in thecontinuum limit and at physical quark masses and we explain our procedure in Section V. In thatsection we define our scaling trajectory, explain how we match the parameters at the differentlattice spacings so that they correspond to the same physics and discuss how we perform the ex-trapolations. We consider this discussion to be a significant component of this paper and believethat this will prove to be a good approach in future efforts to obtain physical results from latticedata. Although we apply the procedures to our data at two values of the lattice spacing, we stressthat the discussion is more general and can be used with data from simulations at an arbitrarynumber of different values of b . In the second half of Section V we then perform the combinedcontinuum–chiral fits in order to obtain our physical results for the decay constants, physical bare01 l +m res (GeV)00.00010.00020.00030.00040.0005 c Q ( G e V ) LOa -1 =1.732 GeVa -1 =2.284 GeV FIG. 54: Topological susceptibility (24 (squares), 32 (circles)). The dashed line is the prediction from LOSU(2) chiral perturbation theory (Eq. (100)) with the chiral condensate computed from the finite volumeLEC’s given in Table XXVII. The solid line denotes the result of the single-parameter fit to the NLOformula given in Eq. (102). quark masses (which are renormalized in Section VI) and for the quantities r and r defined fromthe heavy-quark potential. For the discussion below, it is important to recall that we use the phys-ical pion, kaon and W masses to determine the physical quark masses and the values of the latticespacing and we then make predictions for other physical quantities.In contrast to most other current lattice methods, the DWF formulation gives our simulationsgood control over chiral symmetry, non-perturbative renormalization factors and flavor symmetry.This control allows us to measure and use, as either inputs or predictions: pseudoscalar decayconstants, as well as their ratios; pseudoscalar masses; baryon masses; weak matrix elementsand static potential values, limited only by the statistics achievable for these observables. Theability to predict many observables from the same simulations, provides evidence for the generalreliability of the underlying methods. The good properties of DWF also allow us to test scaling,over this wide range of observables, at unphysical quark masses, since there are no flavor or chiralsymmetry breaking effects to distort a test of scaling. We find scaling violations at the percentlevel, which supports including scaling corrections in only the leading order terms in our light-quark expansions.02As we reduce the quark masses used in the simulations, it is frustrating that there remains a doubtas to the best ansatz to use for the chiral extrapolation. We know of course that for sufficiently light u and d masses the behaviour is given by SU(2) ChPT; what we don’t know is what ”sufficientlylight” means in practice. While in the range of quark masses accessible in our simulations, corre-sponding to 290 - 420 MeV for unitary pions and 225 - 420 MeV for partially quenched pions, ourdata are consistent with NLO SU(2) ChPT, we have seen that they are also consistent with a simpleanalytic ansatz leading to an inherent uncertainty in how best to perform the chiral extrapolation.This is particularly well illustrated in the study of f p , see Fig. 35 for example, where the data iswell represented by all three ans¨atze (including NLO SU(2) ChPT with finite-volume corrections),but the extrapolated values differ as seen in Table XXXI f p = ( ) MeV from the NLO ChPTanalysis with finite-volume corrections and f p = f p = ( )( ) MeV. In Section V E 3 we investigated theincrease in c /dof if the fits are required to pass through the physical value 130.7(4) MeV up tocorrections from lattice artefacts and found c =1.9(7) for the analytic ansatz and an unacceptablylarge value of 5(1) for the NLO ChPT with finite volume corrections. In the future, it will bevery interesting to see how the different ans¨atze for the chiral extrapolation become constrained orinvalidated as we perform simulations with even lighter masses. We point out that the differencein the results from the analyses using the finite-volume ChPT and analytic ans¨atze is much smallerfor the other quantities studied in this paper than for f p .The main physical results of this study are: f p = ( )( ) MeV { Eq. ( ) } ; f K = ( )( ) MeV { Eq. ( ) } ; f K f p = . ( )( ) { Eq. ( ) } ; m MS s ( ) = ( . ± . ) MeV { Eq. ( ) } ; m MS ud ( ) = ( . ± . ) MeV { Eq. ( ) } ; [ S MS ( )] / = ( ) MeV { Eq. ( ) } ; r = . ( ) fm and r = . ( ) fm { Eq. ( ) } . (103)For convenience we also display the equation number where the results were presented earlier inthis paper to help the reader find the corresponding discussion. All the results in Eq. (103) were03obtained after reweighting the strange-quark mass to its physical value at each b , and the renor-malized quark masses were obtained using non-perturbative renormalization with non-exceptionalmomenta as described in Section VI. The low-energy constants obtained by fitting our data toNLO chiral perturbation theory can be found in Sec. V E.The configurations and results presented in this paper are being used in many of our current stud-ies in particle physics phenomenology, including the determination of the B K parameter of neutralkaon mixing in the continuum limit [34]. In parallel to these studies we are exploiting config-urations generated at almost physical pion masses on lattices with a large physical volume ( ∼ D I = / K → pp decay amplitudes were recentlypresented in Refs. [48, 69]. Having access to data with excellent chiral and flavor properties with arange of lattice spacings and quark masses makes this an exciting time indeed for studies in latticephenomenology. Acknowledgments
Appendix A: Separate fits to and data In this section we report on results obtained by fitting the data from the 24 runs at b = .
13 andfrom the 32 runs at b = .
25 separately to the predictions of SU ( ) × SU ( ) ChPT. This comple-ments the material presented in Sections III and IV in which we presented the results for massesand decays constants at each set of quark masses but did not perform the chiral extrapolationsand also that in Section V in which we performed simultaneous chiral and continuum fits to thedata at both lattice spacings. Our main motivation for studying separate fits here is to be ableto compare directly our results obtained with the new data to those in our previous publication[1]. For that reason in this appendix we will be using the same renormalization constant Z A as inour previous publication, which differs from the one used in the global analysis presented in themain part of this paper, see the discussion in Sec. III and App. B for details. We use the samemethod of iterated fits as outlined in our earlier publication [1]; at each lattice spacing we iteratethe combined fits of the meson masses and decay constants with m x ≤ .
01 to the SU ( ) -ChPTformulae, using kaon SU ( ) ChPT to fit the kaon mass and decay constants and the extrapolationin the W -baryon mass until convergence. The pion, kaon, and W masses are used to fix the phys-ical bare quark masses m ud , m s and the lattice scale 1 / a . Predictions for the remaining physicalquantities are then obtained by extrapolation to these physical quark masses. For further detailssee [1]. In the case of the 24 ensembles, the runs have been extended since the publication of[1] (see Sec. II and especially Tab. I for details) so that a direct comparison of the results fromthe previous (smaller) data set with the new extended data set is possible. We quote results fromfits with and without corrections due to finite-volume effects. When including the finite volumecorrections, the terms described in Appendix C of [1] are included in the SU ( ) ChPT in the pionsector (both for the meson masses and decay constants). We also include the correction termscontaining the chiral logarithm of the light quark masses in the kaon decay constant [82] and notethat up to NLO in the light-quark masses, no finite-volume corrections arise in the masses of thekaon and W -baryon. Below we present the physical results in the infinite-volume limit, i.e. after05removing the corrections. Finally, we will perform a na¨ıve continuum extrapolation of the resultsobtained by the separate fits at the two lattice spacings, which can then be compared to resultsfrom the combined chiral-continuum extrapolations using the global fits described in Sec. V. Notethat in this appendix also for the combined chiral-continuum extrapolations we are going to quoteresults obtained using our previous definition of Z A . For that reason the results reported here differslightly from those in the main part of this paper.
1. SU(2)-ChPT fits to data In Tab. XXXIX we summarize our results from the iterative fits to the masses and decay constantsmeasured on the 24 ensembles (see Sec. III for details) and compare them to our earlier results ob-tained with lower statistics [1]. We have performed two kinds of fits: one including the W -baryonmasses determined at all the simulated light-quark masses, m l = . W -baryon masses at the two lightest dynamical quarkmasses m l = .
005 and 0.01 are included. The latter, limited range is also the one used in thecombined chiral-continuum extrapolations in Section V and in the separate fits to the 32 datain the next subsection. In Fig. 55 we plot the combined SU ( ) ChPT fits (without finite-volumecorrections) to the meson masses and decay constants in the pion sector. It is evident that overthe fit range ( m x + m y ) / ≤ .
01, corresponding to a maximum meson mass of about 420 MeV,the data is well described by SU ( ) ChPT. This is also true for the fits including the finite-volumecorrections (not shown).We note that by comparing the results in the first two columns of Tab. XXXIX, which have beenobtained using the same (large) mass-range for the chiral extrapolation of the W -baryon mass, theresults obtained with the increased statistics (for each dynamical light-quark mass the statisticshas nearly been doubled, see Section III) nicely agree with those from our previous publication [1]within the statistical uncertainty. Furthermore, we observe the expected reduction in the statisticalerror. For the remainder of the discussion, we focus on the fits in which only the two lightestdynamical masses have been included in the extrapolation of the W -baryon mass, i.e. the last twocolumns of Tab. XXXIX. The major difference resulting from this change in the fit range is in thevalue of the lattice scale 1 / a , but within 1.4 standard deviations (statistical error only, taking intoaccount correlations) the results still show agreement. Including the finite-volume effects resultsin higher values for the decay constants (both in the chiral limit and at the physical point), which06 Allton et al. [1] increased statisticsno FV-corr. no FV-corr. incl. FV-corr. W : m l ≤ . W : m l ≤ . W : m l ≤ . W : m l ≤ . / a [ GeV ] B MS ( ) [ GeV ] . ( . )( . ) ren . ( . )( . ) ren . ( . )( . ) ren . ( . )( . ) ren f [ MeV ] l l f p [ MeV ] f K [ MeV ] f K / f p m MS ud ( ) [ MeV ] . ( . )( . ) ren . ( . )( . ) ren . ( . )( . ) ren . ( . )( . ) ren m MS s ( ) [ MeV ] . ( . )( . ) ren . ( . )( . ) ren . ( . )( . ) ren . ( . )( . ) ren ˜ m ud : ˜ m s aB a f L ( ) × L ( ) × ( L ( ) − L ( ) ) × -0.71(0.62) -0.09(0.45) 0.03(0.45) 0.01(.49) ( L ( ) − L ( ) ) × a ˜ m ud a ˜ m s ( ) ChPT fits to the 24 data (without and with finite-volume correc-tions) compared to those from [1] obtained with lower statistics (without finite-volume corrections). We alsoquote in the lower part of the table the SU ( ) ChPT fit parameters aB , a f , L ( ) i (at the scale L c = a ˜ m ud , s in lattice units. Only statistical uncertainties are quoted except for quarkmasses and the LEC B renormalized in the MS-scheme at 2 GeV where also the systematic uncertaintyfrom the renormalization constant is quoted. (Mass renormalization constant at 1 / a = . ( ) GeV: Z m = . ( . ) stat ( . ) ren and at 1 / a = . ( ) GeV: Z m = . ( . ) stat ( . ) ren .) f xy m y combined f p , m p fitm l = 0.005 m x =0.001m x =0.005m x =0.010m x =0.020m x =0.030m x =0.040 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 0 0.01 0.02 0.03 0.04 m xy / ( m a v g + m r e s ) m y combined f p , m p fitm l =0.005 m x =0.001m x =0.005m x =0.010m x =0.020m x =0.030m x =0.040 0.080 0.085 0.090 0.095 0.100 0.105 0.110 0 0.01 0.02 0.03 0.04 f xy m y combined f p , m p fitm l = 0.01 m x =0.001m x =0.005m x =0.010m x =0.020m x =0.030m x =0.040 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 0 0.01 0.02 0.03 0.04 m xy / ( m a v g + m r e s ) m y combined f p , m p fitm l =0.01 m x =0.001m x =0.005m x =0.010m x =0.020m x =0.030m x =0.040 FIG. 55: Combined SU ( ) ChPT fits (without finite-volume corrections) for the meson decay constants (leftcolumn) and masses (right column) on the 24 data set at m l = . (top row) and 0.01 (bottom row) . Onlypoints marked with circles , corresponding to the range ( m x + m y ) / ≤ .
01 are included in the fits. is a statistically significant effects (taking the correlations into account). In Tab. XL we comparethe decay constants and their ratio obtained from the separate fits with the corresponding resultsfrom the global analysis at the simulated, finite value of the lattice spacing (i.e. not extrapolatedto the continuum, see Sec. V and especially Tabs. XXXI, XXXII, XXXIII but note the differencedue to the use of our previous definition of Z A here). We are reassured by the observed agreementbetween the results obtained using the global fits with those obtained using our previous strategy08 f p [ MeV ] f K [ MeV ] f K / f p no FV-corr. 24 , b = .
13 separate 124.4(3.6) 151.0(3.7) 1.214(0.012)global 123(2) 150(2) 1.215(0.009)32 , b = .
25 separate 120.4(1.9) 147.1(2.0) 1.222(0.007)global 121(2) 147(2) 1.222(0.006)incl. FV-corr. 24 , b = .
13 separate 126.4(3.6) 152.1(3.7) 1.204(0.012)global 126(2) 151(2) 1.204(0.009)32 , b = .
25 separate 122.3(1.9) 148.1(2.0) 1.212(0.007)global 123(2) 149(2) 1.210(0.006)TABLE XL: Comparison of the pion and kaon decay constants and their ratios at finite lattice spacing fromseparate (see Tabs. XXXIX, XLI) and global fits using our previous definition of Z A . in Ref. [1] which was developed at that time to describe data at only a single lattice spacing.
2. SU(2)-ChPT fits to data The results of a separate fit on the 32 data set are summarized in Tab. XLI. Here we only includedthe W -baryon masses from the m l = . ( m x + m y ) / ≤ . ( ) ChPT.As was already the case for the 24 ensembles, taking finite-volume corrections into account alsoleads to a good description of the data and results in higher values for the decay constants at thephysical point and in the chiral limit. Again, taking the correlations into account, we note that thisis a statistically significant effect. As was also the case on the 24 ensembles, we observe a goodagreement for the decay constants and their ratio between the results of the separate fits to the 32 data and the results from the global fits at finite lattice spacing, see Tab. XL.09 no FV-corr. FV-corr. incl.1 / a [ GeV ] B MS ( ) [ GeV ] . ( . )( . ) ren . ( . )( . ) ren f [ MeV ] l l f p [ MeV ] f K [ MeV ] f K / f p m MS ud ( ) [ MeV ] . ( . )( . ) ren . ( . )( . ) ren m MS s ( ) [ MeV ] . ( . )( . ) ren . ( . )( . ) ren ˜ m ud : ˜ m s aB a f L ( ) × -0.75(0.79) -1.21(.82) L ( ) × ( L ( ) − L ( ) ) × -0.93(0.42) -1.03(0.45) ( L ( ) − L ( ) ) × a ˜ m ud a ˜ m s ( ) ChPT fits to the 32 data (without and with finite-volume corrections).We also quote in the lower part of the table the SU ( ) ChPT fit parameters aB , a f , L ( ) i (at the scale L c = a ˜ m ud , s in lattice units. Only statistical uncertainties are quoted except forquark masses and the LEC B renormalized in the MS-scheme at 2 GeV where also the systematic uncertaintyfrom the renormalization constant is quoted. (Mass renormalization constant at 1 / a = . ( ) GeV: Z m = . ( . ) stat ( . ) ren .) f xy m y combined f p , m p fitm l = 0.004 m x =0.002m x =0.004m x =0.006m x =0.008m x =0.025m x =0.030 3.3 3.4 3.5 3.6 3.7 3.8 3.9 0 0.01 0.02 0.03 m xy / ( m a v g + m r e s ) m y combined f p , m p fitm l =0.004 m x =0.002m x =0.004m x =0.006m x =0.008m x =0.025m x =0.030 0.060 0.065 0.070 0.075 0.080 0 0.01 0.02 0.03 f xy m y combined f p , m p fitm l = 0.006 m x =0.002m x =0.004m x =0.006m x =0.008m x =0.025m x =0.030 3.3 3.4 3.5 3.6 3.7 3.8 3.9 0 0.01 0.02 0.03 m xy / ( m a v g + m r e s ) m y combined f p , m p fitm l =0.006 m x =0.002m x =0.004m x =0.006m x =0.008m x =0.025m x =0.030 0.060 0.065 0.070 0.075 0.080 0 0.01 0.02 0.03 f xy m y combined f p , m p fitm l = 0.008 m x =0.002m x =0.004m x =0.006m x =0.008m x =0.025m x =0.030 3.3 3.4 3.5 3.6 3.7 3.8 3.9 0 0.01 0.02 0.03 m xy / ( m a v g + m r e s ) m y combined f p , m p fitm l =0.008 m x =0.002m x =0.004m x =0.006m x =0.008m x =0.025m x =0.030 FIG. 56: Combined SU ( ) ChPT fits (without finite-volume corrections) for the meson decay constants (leftcolumn) and masses (right column) on the 32 data set at m l = . (top row) , 0.006 (middle row) , and0.008 (bottom row) . Only points marked with circles , corresponding to the range ( m x + m y ) / ≤ .
008 areincluded in the fits.
3. Extrapolation to the Continuum Limit
With the results obtained from separate chiral extrapolations on the 24 (extended statistics) andthe 32 data sets (see the two previous subsections, respectively) we can perform a na¨ıve contin-uum limit extrapolation assuming a -scaling. Of course, with only two lattice spacings available,we are not able to confirm this scaling behaviour. Further caveats include the fact that here, forsimplicity, we did not use reweighting and so the dynamical strange-quark mass is not tuned toexactly the same value on the two data sets and indeed is not exactly the physical one on eitherset. Also, the dynamical light-quark mass ranges are a little different at the two lattice spac-ings, corresponding to unitary pion masses in the range 330–420 MeV on the coarser 24 latticesand 290–400 MeV on the finer 32 lattices (a similar statement is true for the partially-quenchedmasses). One might therefore expect a larger uncertainty in the chiral extrapolation of the 24 results. In the na¨ıve continuum ansatz followed here, we are not taking into account this effect.Because of this, and maybe more importantly, since two separate chiral extrapolations have beenperformed (one at each of the two values of the lattice spacing), the continuum extrapolation isnot completely disentangled from the chiral extrapolation. Recall that in our procedure for theglobal fits described in the main part of this paper, these two extrapolations are indeed disentan-gled. There this is achieved by adding O ( a ) terms into the two functions, such that the chiral andcontinuum extrapolations are performed simultaneously and independently from each other.In Tab. XLII we repeat the results obtained at the two different lattice spacings (with and with-out finite-volume corrections) and give the values extrapolated to the continuum limit assuming a scaling. Fig. 57 illustrates the continuum extrapolation of the various quantities (only resultsobtained without taking into account finite-volume corrections are shown there). Note, that thetwo points at the different lattice spacings are completely uncorrelated, the only correlation in thedata for the continuum extrapolation is between the uncertainty in the lattice spacing (the “x”-datum) and the quantity itself at that lattice spacing (the “y”-datum). These correlations weretreated by the super-jackknife method which we have been using in our earlier work and which isclearly explained in [73, 74]. For comparison, Tab. XLII also contains our results from the com-bined continuum-chiral extrapolation as described in the main part of this paper but here using ourprevious definition of Z A . As one can see, the combined continuum-chiral extrapolation gives asubstantially smaller (up to a factor of 5) statistical uncertainty compared to the na¨ıve continuumextrapolation. The main reason, of course, is the correlation in the combined fits between the two12 no FV-corr.separate fits na¨ıve CL comb. chiral/CL24 , b = .
13 32 , b = . a [ fm ] → → f [ MeV ] l l f p [ MeV ] f K [ MeV ] f K / f p , b = .
13 32 , b = . a [ fm ] → → f [ MeV ] l l f p [ MeV ] f K [ MeV ] f K / f p and 32 data sets ( W masses from m l ≤ . data set, cf. Tabs. XXXIX and XLI) and their na¨ıve continuum limit assuming a -scaling (see Fig. 57)compared to results from the combined chiral-continuum extrapolation using the previous definition of Z A .The top table contains results without finite-volume corrections whereas the results in the bottom table wereobtained by including finite-volume effects.
100 105 110 115 120 125 130 0 0.089 a [fm ]f p [MeV] 130 135 140 145 150 155 0 0.089 a [fm ]f K [MeV] 90 95 100 105 110 115 120 0 0.089 a [fm ]f [MeV] 1.2 1.21 1.22 1.23 1.24 1.25 1.26 1.27 0 0.089 a [fm ]f K /f p a [fm ]bar l a [fm ]bar l FIG. 57: Results from separate fits (without finite-volume corrections) to the 24 and 32 data sets (blackpoints) and the na¨ıve continuum-limit extrapolation (blue asterisks) for selected quantities assuming a -scaling. For details see Subsec. A 3 and Tab. XLII. O ( a ) corrections for the leading-orderterms, as is consistent with our power counting scheme. In this way, the continuum extrapo-lation in the combined fits is also more constrained, leading to a smaller statistical uncertainty.Comparing the results of the na¨ıve continuum extrapolation and the combined continuum-chiralextrapolation for the quantities in Tab. XLII we observe agreement better than 0.5- s (taking intoaccount correlations) for all quantities except for ¯ l , where the agreement still holds at the 1- or1.5- s level (without and with taking FV-corrections into account, respectively). It is reassuring,that the results from the two methods agree well, although the value of this statement is limited,given the large (statistical) uncertainty of almost 10% for the decay constants or even more incase of the LECs from the na¨ıve method. However, it should be noted that the same agreementholds, not only for the continuum values, but also for the results obtained in the separate fits ascompared to the predictions of the global fit made for the finite lattice spacings. This has alreadybeen discussed in the previous subsections and is shown in Tab. XL. Appendix B: Determining Z A As pointed out by Sharpe [17] and refined in Ref. [1], the normalization of the partially conservedaxial current defined for domain wall fermions [75] is expected to deviate from that of the con-ventionally normalized continuum current by an amount of order m res a . Here and below whenmaking such estimates we will introduce the explicit lattice spacing a and express the residualmass in physical units in order to make the comparison of various terms in a Symanzik expan-sion in powers of a easier to recognize. Since such a deviation can be viewed as O ( ma ) whichis formally larger than the O ( ma ) which we neglect in our power counting scheme and becausethe normalization of this axial current plays a central role in our determination of the importantquantities f p and f K , we have calculated this normalization factor Z A numerically. We explain ourmethod and result in this appendix. The first subsection contains a discussion of the theoreticalissues and explains the basis for our method of determining Z A . The second subsection describesthe actual calculation and results.15
1. Determining the normalization of A m To determine the normalization of A m we compare the matrix element of four distinct domain wallfermion currents. The first two are the conserved/partially conserved vector and axial currents V a m ( x ) and A a m ( x ) respectively, where a and m are flavor and space-time indices. These currentswere introduced by Furman and Shamir [75] and involve fermion fields evaluated on each of the L s x and x + ˆ e m where ˆ e m is a unit vectorpointing the m th direction. Thus, these currents are local but distributed in the fifth dimension andone-link non-local in space-time. While this vector current is exactly conserved, the divergence ofthe axial current contains the usual mass term and a mid-point term J a q . In the long-distance limitthis midpoint term can be decomposed into the residual mass term, a piece that is convenientlywritten as ( − Z A ) times the divergence of the same axial current and a final term of dimensionfive which we write out explicitly as the sum of the dimension-five, chiral rotation of the usualclover term and the four-dimensional Laplacian applied to the pseudoscalar density: J a q = m res q g l a q + − Z A D m A a m + c q s mn F mn l a q + c ¶ m ¶ m q g l a q . (B1)In Equation (B1) l a is the generator which acts on the fermion fields corresponding to the flavorindex a while q ( x ) and q ( x ) are the “physical”, four-dimensional quark fields obtained by evaluat-ing the five-dimensional domain wall fields on the s = s = L s − V a m ( x ) and A a m ( x ) , constructed in the standard way from the four-dimensional quark fields, q ( x ) and q ( x ) . These currents are localized in all five dimensions and neither is conserved.Finally it will also be convenient to introduce the scalar densities q ( x ) q ( x ) , q ( x ) l a q ( x ) from whichthe domain fermion mass is constructed and their chiral transforms q ( x ) g q ( x ) , q ( x ) l a g q ( x ) .These four classes of operators will be labeled S ( x ) , S a ( x ) , P ( x ) and P a ( x ) .Following Symanzik, we can add improvement terms to each of these six operators to insure thattheir Green’s functions, when evaluated with an appropriately improved action, will agree withthe corresponding continuum Green’s functions up to errors of order a n . For our present purposes,accuracy up to O ( am ) where m is a quark mass in physical units, will be sufficient. Since m res and m have a similar size, we are explicitly attempting to control the m res a corrections describedabove. We do not attempt to explicitly remove O ( a ) terms since these will be eliminated by thefinal linear extrapolation a → V S a m , A S a m , S S a and P S a for the Symanzik-improved vector current, axial cur-rent, scalar density and pseudoscalar density respectively, keeping improvement terms which arenominally of order a and imposing charge conjugation symmetry, we find: V S a m = Z V V a m + C V ¶ n q s mn l a q (B2) A S a m = Z A A a m + C A ¶ m P a (B3) V S a m = Z V V a m + C V ¶ n q s mn l a q (B4) A S a m = Z A A a m + C A ¶ m P a (B5) S S a = Z S S a (B6) P S a = Z P P a . (B7)In contrast to the Symanzik-improved current operators, we have not specified a normalizationconvention for the operators S S a and P S a . Adopting definitive conventions for S S a and P S a is notneeded here beyond the requirement that those conventions are consistent with S S a ± P S a belongingto the ( , ) / ( , ) representations of the SU(3) L × SU(3) R flavor symmetry.Because the operators S and P contain no vector indices, any correction terms must increase thedimension by two and we have chosen to neglect such O ( a ) contributions. Thus, Eqs. (B6) and(B7) are particularly simple. However, we can also drop the dimension four, O ( a ) correction termsto Eqs. (B2)-(B5). This can be established by considering the chiral structure of the Symanzikand conserved/partially conserved current operators. Ignoring effects of order m , the Symanzikcurrents will couple to pairs of quarks which are either left- or right-handed. Likewise the domainwall conserved/partially conserved current operators couple to a pair of quarks with the samevalue of the coordinate s in the fifth dimension. For s = s = L s − s moves into the fifth-dimensional bulk, the17amplitude for coupling to such physical modes decreases until when s ≈ L s / m res a . Of course, the s ≈ s ≈ L s − s = L s −
1. Since the four, dimension-four improvement terms included in Eqs. (B2)-(B5) involvepairs of quarks with opposite handedness, such terms require a complete propagation across thefifth dimension if they are to couple to the conserved/partially conserved or local currents. This istrue even for the terms with general s which appear in the former currents. Thus, these correctionterms involve an additional power of m res a and are of order m res a and can be neglected in ourpower counting scheme.With this simplification, we can demonstrate that to this order the following relations hold: Z V = Z V = Z A (B9) Z S = Z P . (B10)Equation (B8) follows easily from the fact that V a m is conserved at finite lattice spacing and hasbeen given the conventional normalization. Equations (B9) and (B10) can each be shown usingessentially the same argument which we will now review.In the massless continuum theory the operators q c l a g m ( ± g ) q c are independent involving onlyright-handed/left-handed degrees of freedom. Here the label c indicates continuum . This impliesthe vanishing of the Symanzik-improved Green’s function: D ( V S a m + A S a m )( x )( V S a n − A S a n )( y ) E = . (B11)This same property is obeyed by the local domain wall currents up to order ( m res a ) since non-vanishing terms which can contribute to the DWF version of Eq. (B11) must connect both fermiondegrees of freedom between the left and right walls requiring two-traversals of the fifth dimensionand hence are of order ( m res a ) [17, 76]. It is then easy to see that these two behaviors can beconsistent through order m res a only if Z V = Z A through order m res a . We need only examine the18mixing between V S a m ± A S a m that is generated by Z V − Z A : D ( V S a m + A S a m )( x ) · ( V S a n − A S a n )( y ) E (B12) = D ( Z V V a m + Z A A a m )( x ) · ( Z V V a n − Z A A a n )( y ) E = D h ( Z V + Z A )( V a m + A a m )( x ) + ( Z V − Z A )( V a m − A a m )( x ) i · h ( Z V + Z A )( V a n − A a n )( y ) + ( Z V − Z A )( V a n + A a n )( y ) i E . The product of the left-most operators in the square brackets on the right-hand side of Eq. (B12)cannot mix at order m res because of their construction from domain wall quark fields as explainedabove. Likewise the product of the right-most terms also vanishes. However, the two cross termshave non-zero correlators implying that for the entire expression to be of order m , the difference Z V − Z A must be of order ( m res a ) , demonstrating the intended result. A very similar argumentcan be constructed which shows that Z S = Z P through order m res a . One must invoke the flavorstructure and, for example, consider correlators between ( S − iS )( x ) + ( P − iP )( x )) and ( S + iS )( y ) + ( P + iP )( y )) which also must vanish in the chiral limit. Here a = , a = − Z A .Since at low energies the left- and right-hand sides of Eqs. (B4) and (B5) must have identicalmatrix elements, the ratio of long-distance correlators computed with the Symanzik and localcurrents must give identical constants: Z V = Z A . Thus, we have established: (cid:10) V S ai ( x ) V ai ( y ) (cid:11)(cid:10) V ai ( x ) V ai ( y ) (cid:11) = (cid:10) A S a ( x ) P a ( y ) (cid:11)(cid:10) A a ( x ) P a ( y ) (cid:11) (B13)where we have introduced the fixed spatial index i , the temporal index 0 and sources V ai ( y ) and P a ( y ) that will correspond to those used in our actual calculation. Next we can use the long-19distance equality represented by Eqs. (B2) and (B3) to write1 = (cid:10) V S ai ( x ) V ai ( y ) (cid:11)(cid:10) V ai ( x ) V ai ( y ) (cid:11) (B14) Z A = (cid:10) A S a ( x ) P a ( y ) (cid:11)(cid:10) A a ( x ) P a ( y ) (cid:11) . (B15)Then we can combine Eqs. (B13), (B14) and (B15) to yield an equation for Z A which does notinvolve the Symanzik currents: Z A = (cid:10) A a ( x ) P a ( y ) (cid:11)(cid:10) A a ( x ) P a ( y ) (cid:11) · h V ai ( x ) V ai ( y ) i (cid:10) V ai ( x ) V ai ( y ) (cid:11) , (B16)which determines Z A in terms of four correlators which we have evaluated directly in our latticecalculation.In order to relate the discussion of the Symanzik improved operators given in Eqs. (B2)-(B7) withthe operators appearing in Eq. (B1), we should recognize that the quantity Z A has been introducedin two places. The most important is in the relation between the Symanzik current and the partiallyconserved domain wall operator in Eq. (B3). It is this quantity that is determined in Eq. (B16) andwhich is needed to give a physical normalization to the axial current matrix elements determined inour calculation. However, the quantity Z A also appears in the expression for J q given in Eq. (B1).For completeness, we will now demonstrate that these two quantities are in fact the same up toorder ( m res a ) .This is easily done by introducing a flavor-breaking mass term qMq into the DWF action, exam-ining the divergence equations obeyed by V a m and A a m and using the relation Z S = Z P establishedabove. With the additional mass term the conserved/partially conserved vector and axial currentsobey the lattice divergence equations, through O ( m res a ) : D m V a m = q [ l a , M ] q (B17) D m A a m = q { l a , M } g q + m res q g q − ( Z A − ) D m A a m . (B18)Taking the Z A − S a and P a are symmetrically normalized ( Z S = Z P ) , we can conclude that the operators V a m and Z A A a m must be related to the corresponding Symanzik currents by the same factor. Thisestablishes that our two definitions of Z A are consistent.We will conclude this analysis with a brief discussion of the effects of the explicit quark mass, m f , on the operator product expansion represented by Eq. (B1) and on the Symanzik-improved20operators given in Eqs. (B2)-(B7). Although m f explicitly connects the s = s = L s − J q appearing on the left hand side of Eq. (B1)to create effects with arbitrary chiral properties. Thus, we expect multiplicative corrections ofthe form ( + b i m f a ) ≤ i ≤ to each of the four terms on the right hand side of Eq. (B1). In thecase of the left-most term the correction is of order m f m res a while for the remaining three termsthe corrections are of order m f m res a or m f m res a , all beyond the level of accuracy of the currentpaper. The conclusion that Z V = m res a (and order m f a ) prevents the appearanceof a factor 1 + b ( m f a ) multiplying the Z V in Eq. (B2). The argument that Z A = Z V and Z S = S P with corrections of order ( m res a ) applies equally well to the left-right mixings created by m f butagain the allowed m f m res a and ( m f a ) terms are negligible within our present power countingscheme so Eqs. (B4)-(B7) need no O ( m f a ) corrections. Lastly, consider adding a factor of theform ( + b ( m f a )) multiplying the Z A on the right-hand side of Eq. B3. As explained above, asimilar correction to Z A appearing in Eq. (B1) carries the additional suppression of one powerof m res a . Since the equality derived above between the Z A factors appearing in the divergenceequation, Eq. (B1), and the Symanzik-improved current A a m , in Eq. (B3), holds at order m f a sucha 1 + b ( m f a ) factor is not allowed in Eq. (B3). Thus, no m f a terms need to be introduced into theequations presented in this appendix.
2. Computational method and results
We have evaluated the two factors in Eq. (B16) to determine Z A on both the 32 × b = . m l = . × b = .
13 ( m l = . Z A / Z A duplicate those from the calculation of Z A described in Sections IIIand IV. In this appendix we add the factor Z A in the denominator because we are now determiningthe deviation of this factor from unity. We do not simply use the results presented earlier in thepaper because our calculation of Z V / Z V has been performed on a subset of the configurationsanalyzed earlier and results for Z A / Z A are needed on this same subset of configurations if ratioswith meaningful jackknife errors are to be determined.The ratio Z A / Z A was computed from the same ratio of current-pseudoscalar correlators studied inSections III and IV, using the method specified in Ref. [77]. Similar methods are used to compute21 b m l Z A / Z A Z V / Z V Z V / Z A Fit range N meas − m res − m res Z A / Z A , Z V / Z V and Z V / Z A computed on six ensembles. The rowswith quark mass − m res contain the chiral extrapolation to the light quark mass m l = − m res . The left-handportion of the fit range gives that used for the axial current ratio while the right hand portion that for thevector current. For the Z V / Z V calculation the data at t and 63 − t were combined for 0 ≤ t < Z V / Z V using the ratio of vector correlators Z V Z V = (cid:229) i = (cid:229) ~ x D V ai ( ~ x , t ) V ai ( ~ , ) E (cid:229) i = (cid:229) ~ x D V ai ( ~ x , t ) V ai ( ~ , ) E , (B19)an equation expected to be valid for time separations t much larger than one lattice spacing: t ≫ a .Figure 58 shows the right-hand side of Eq. (B19) as a function of time for the case of the lightestmass for each of the 32 and 24 ensembles. A constant fit to plateau regions identified by thehorizontal lines was then used to determine the Z V / Z V on the left-hand side of this equation.Fig. 59 displays the chiral extrapolation of the two quantities Z A / Z A and Z V / Z V on both sets ofensembles.Two useful results follow from this Appendix. First the ratio Z V / Z A differs from unity on ourtwo ensembles and that difference decreases more rapidly than a with increasing b . Thus, wewill obtain more accurate results in our continuum extrapolation from both matrix elements ofthe local axial current and our NPR calculations which are normalized using off-shell Green’sfunctions containing the local vector and axial currents if we convert the normalization of theselocal currents to the usual continuum normalization by using the ratio Z V / Z V instead of the ratio Z A / Z A , the quantity which we have used in previous work for such conversions. The values of Z V / Z V presented in Table XLIII are therefore used to normalize the results presented in the current22 FIG. 58: Plots of the correlator ratio which determines the renormalization factor Z V / Z V as a function oftime. The left panel shows results from the 32 , m l = .
004 ensemble while the right panel the result fromthe 24 , m l = .
005 ensemble. The horizontal line with error bands in each panel shows the fitting rangeand the result obtained in each case.FIG. 59: The quantities Z A / Z A and Z V / Z V extrapolated to the chiral limit for the 32 (left panel) and 24 (right panel) ensembles. paper and are the second result obtained in this appendix. Because these ratios were calculatedon a smaller subset of configurations than were used for our main results, we have included theirstatistical fluctuations as independent within our superjackknife, statistical error analysis. Sincethese fluctuations are at or below the 0.5% level, this omission of possible statistical correlationsis unimportant.23 Appendix C: Statistical errors of reweighted quantities
In this appendix we discuss the statistical errors that should be expected when Monte Carlo datais reweighted to obtain results for a gauge or fermion action that is different from that used togenerate the data. Throughout this discussion we will make the assumption that the reweightingfactors are not correlated with the data. Of course, if this assumption were exactly true then thereweighting would not be needed. However, the correlation between the data and reweightingfactors is often small in practice and neglecting this correlation may well provide a reasonablyaccurate view of the resulting errors. As we will show, with this assumption the usual analysis ofthe statistical errors applies easily to reweighted data and yields simple, useful formula which wepresent here.Consider a quantity x and the corresponding ordered ensemble of N Monte Carlo configurationswith corresponding measured values { x n } , 1 ≤ n ≤ N . For each of these N configurations we willdetermine a reweighting factor w n so that the final, reweighted quantity of interest is given by h x i N = (cid:229) Nn = x n w n (cid:229) Nn = w n . (C1)Here the single brackets h . . . i N indicate an average over a single Monte Carlo ensemble of N samples. In this appendix we are interested in how the statistical fluctuations in the quantity h x i N are affected by the operation of reweighting. We can then express the true value for x N as x N = DD h x i N EE (C2)where the double brackets hh . . . ii indicate a “meta” average over many equivalent Monte Carloensembles. The statistical fluctuation present in a particular result h x i N can then be characterizedby the average fluctuation of h x i N about x N :Error ( x ) = rDD ( h x i N − x N ) EE . (C3)A quantity such as h x N i , defined in Eq. (C1) as a ratio of averages, will be a biased estimatorof the physical result which must be determined in the limit N → ¥ . Thus, the meta average x N = hhh x i N ii will differ from the true result by terms of order 1 / N . While these 1 / N correctionsare not difficult to enumerate and estimate from our data, these corrections are not the subject ofthe present appendix and will not be considered further here. Instead we will study how the sizeof the statistical fluctuations of h x N i about x N is affected by the reweighting. Thus, the quantity24Error ( x ) defined in Eq. (C3) describes the average deviation of h x N i from x N not from the N → ¥ limit of x N .We will now work out an expression for Error ( x ) in the case that nearby measurements x n and x n + l in a single Markov chain (or reweighting factors w n and w n + l ) are correlated but with theassumption that x n and w n + l are not: DD(cid:16) h x i N − x N (cid:17) EE = **(cid:18) (cid:229) Nn = x n w n (cid:229) n w n − x N (cid:19) (cid:229) Nn ′ = x n ′ w n ′ (cid:229) n ′ w n ′ − x N !++ (C4) = ** (cid:0) (cid:229) Nn = x n w n − x N (cid:229) Nn = w n (cid:1) (cid:0) (cid:229) Nn ′ = x n ′ w n ′ − x N (cid:229) Nn ′ = w n ′ (cid:1)(cid:0) (cid:229) Nn = w n (cid:1) (cid:0) (cid:229) Nn ′ = w n ′ (cid:1) ++ (C5) = ** (cid:0) (cid:229) Nn = ( x n − x N ) w n ) (cid:1) (cid:0) (cid:229) Nn ′ = ( x n ′ − x N ) w n ′ (cid:1)(cid:0) (cid:229) Nn = w n (cid:1) (cid:0) (cid:229) Nn ′ = w n ′ (cid:1) ++ (C6) = (cid:229) Nn = (cid:229) N − nl = − n nDD ( x n − x N )( x n + l − x N ) EEDD w n w n + l EEoDD (cid:229) Nn = w n EE , (C7)where in the last line we have used our assumption of the lack of correlation between the x n and w n to write the average of their product as the product of their separate averages. We have alsoassumed that our sample size N is sufficiently large that correlated fluctuations of the averages inthe numerator and denominator will be sufficiently small that the average of the original ratios andproducts can be replaced by the corresponding ratios and products of the individual averages.This result can be cast in a simple form if we define the three averages: d x = DD ( x n − x N ) EE (C8) w = DD w n EE (C9) w = DD w n EE (C10)(where d x is the usual width of the distribution of the measured quantity x n ) and the two autocor-relation functions: C ( l ) = DD ( x n − x N )( x n + l − x N ) EE d x (C11) W ( l ) = DD w n w n + l EE w , (C12)defined so that C ( ) = W ( ) =
1. Making the conventional assumption that the range of l overwhich the correlation function C ( l ) is non-zero is small compared to the sample size N and using25the quantities defined above, we can rewrite Eq. (C7) as DD(cid:16) h x i N − x N (cid:17) EE = d x (cid:229) L max l = − L max C ( l ) W ( l ) w N ( w ) (C13) = d x t corr N w ( w ) (C14)where the autocorrelation time t corr is defined as t corr = L max (cid:229) l = − L max C ( l ) W ( l ) . (C15)The limit L max is chosen to be larger than the region within which C ( l ) is non-zero and has beenintroduced as a reminder that when working with a single finite sample, one must take care toevaluate the limit of large N before the limit of large L max . Finally, Eq. (C14) can be written in theconventional form Error ( x ) = s d x N eff (C16)where the effective number of configurations N eff is given by: N eff = N t corr w w . (C17)This result makes precise a number of aspects of reweighting that may be useful to understand.In the case that there are no autocorrelations so t corr =
1, the ratio w / w expresses the degree towhich the reweighting process selectively samples the original data and degrades the initial statis-tics. The general inequality w / w ≤ w n do not vary with n . In the extreme case that asingle sample w n dominates the averages then w / w = / N and N eff =
1. Thus, in the case ofuncorrelated data (which is the case for most of the results presented here) we should expect thestatistical fluctuations to grow as the degree of reweighting increases by the factor w / w .Including autocorrelations makes the effects of reweighting on the size of the statistical fluctu-ations less certain because the behavior of the factors 1 / t corr and w / w in Eq. (C17) becomeentangled. In the limit in which the autocorrelation time associated with the measured quantity x n alone, t x = L max (cid:229) l = − L max C ( l ) , (C18)26becomes much larger than that of the reweighting factor w n , then the majority of the sum inEq. (C15) contributing to t corr will come from values of l where DD w n w n + l EE ≈ DD w EE so that t corr ≈ t x w w . (C19)In this case the error given by Eq. (C16) reduces to the standard expression p d x t x / N that holdsif no reweighting is performed! Of course, this is easy to understand. When such long autocorre-lation times are involved, the average over the autocorrelation time is providing an average overthe reweighting factors w n which is sufficiently precise that the error-enhancing fluctuations inthe reweighting factors are averaged away. Given the large size of the fluctuations between thereweighting factors and the relatively short autocorrelation times seen in our data, it is unlikelythat this averaging would be seen in the results presented here.A second type of behavior for t corr occurs if the w n are relatively uncorrelated and w ≫ w so thatonly the l = t corr =
1. In this case reweightinghas removed the effects of autocorrelation but increased the statistical fluctuations by the factor w / w which was assumed to be large. Here the fluctuation-enhancing effects of autocorrelationsand reweighting are not compounded. [1] C. Allton et al. (RBC-UKQCD), Phys. Rev. D78 , 114509 (2008), 0804.0473.[2] D. J. Antonio et al. (RBC and UKQCD), Phys. Rev.
D75 , 114501 (2007), hep-lat/0612005.[3] C. Allton et al. (RBC and UKQCD), Phys. Rev.
D76 , 014504 (2007), hep-lat/0701013.[4] D. J. Antonio et al. (RBC), Phys. Rev. Lett. , 032001 (2008), hep-ph/0702042.[5] P. A. Boyle et al. , Phys. Rev. Lett. , 141601 (2008), 0710.5136.[6] P. Boyle, J. Flynn, A. Juttner, C. Kelly, C. Maynard, et al. , Eur. Phys. J. C (2010), arXiv:1004.0886.[7] T. Yamazaki et al. (RBC+UKQCD), Phys. Rev. Lett. , 171602 (2008), 0801.4016.[8] T. Yamazaki et al. , Phys. Rev.
D79 , 114505 (2009), 0904.2039.[9] Y. Aoki et al. , Phys. Rev.
D82 , 014501 (2010), 1003.3387.[10] Y. Aoki et al. (RBC-UKQCD), Phys. Rev.
D78 , 054505 (2008), 0806.1031.[11] N. Christ, C. Dawson, T. Izubuchi, C. Jung, Q. Liu, et al. (2010), * Temporary entry *,arXiv:1002.2999.[12] C. Albertus et al. , Phys. Rev.
D82 , 014505 (2010), 1001.2023. [13] Y. Aoki et al. , Phys. Rev. D78 , 054510 (2008), 0712.1061.[14] C. Sturm et al. , Phys. Rev.
D80 , 014501 (2009), 0901.2599.[15] M. Gorbahn and S. Jager (2010), 1004.3997.[16] L. G. Almeida and C. Sturm, Phys.Rev.
D82 , 054017 (2010), arXiv:1004.4613.[17] S. R. Sharpe (2007), 0706.0218.[18] D. J. Antonio et al. (RBC and UKQCD) (2007), arXiv:0705.2340 [hep-lat].[19] Y. Shamir, Nucl. Phys.
B406 , 90 (1993), hep-lat/9303005.[20] Y. Iwasaki UTHEP-118.[21] M. Hasenbusch and K. Jansen, Nucl. Phys.
B659 , 299 (2003), hep-lat/0211042.[22] T. Takaishi and P. de Forcrand, Phys. Rev.
E73 , 036706 (2006), hep-lat/0505020.[23] C. Urbach, K. Jansen, A. Shindler, and U. Wenger, Comput. Phys. Commun. , 87 (2006), hep-lat/0506011.[24] P. A. Boyle, Computer Physics Communications , 2739 (2009).[25] Y. Aoki et al. , Phys. Rev.
D72 , 114505 (2005), hep-lat/0411006.[26] M. Luscher and F. Palombi, PoS
LATTICE2008 , 049 (2008), 0810.0946.[27] A. Hasenfratz, R. Hoffmann, and S. Schaefer, Phys. Rev.
D78 , 014515 (2008), 0805.2369.[28] C. Jung (2010), 1001.0941.[29] H. Ohki et al. (2009), 0910.3271.[30] S. Aoki et al. (PACS-CS), Phys. Rev.
D81 , 074503 (2010), 0911.2561.[31] R. Baron et al. (ETM), PoS
LATTICE2008 , 094 (2008), 0810.3807.[32] T. Ishikawa, Y. Aoki, and T. Izubuchi, PoS
LAT2009 , 035 (2009), 1003.2182.[33] K. Ogawa and S. Hashimoto, Prog. Theor. Phys. , 609 (2005), hep-lat/0505017.[34] P. Boyle et al. (RBC-UKQCD Collaboration) In Preparation.[35] T. Blum et al. (RBC), Phys. Rev.
D68 , 114506 (2003), hep-lat/0110075.[36] E. Eichten, K. Gottfried, T. Kinoshita, K. D. Lane, and T.-M. Yan, Phys. Rev.
D17 , 3090 (1978).[37] C. Alexandrou et al. , Nucl. Phys.
B414 , 815 (1994), hep-lat/9211042.[38] F. Berruto, T. Blum, K. Orginos, and A. Soni, Phys. Rev.
D73 , 054509 (2006), hep-lat/0512004.[39] T. Blum (RBC and UKQCD), PoS
LATTICE2008 , 096 (2008).[40] C. M. Maynard (RBC and UKQCD), PoS
LAT2009 , 091 (2009), 1001.5203.[41] S. N. Syritsyn et al. , Phys. Rev.
D81 , 034507 (2010), 0907.4194.[42] L. Lellouch, PoS
LATTICE2008 , 015 (2009), 0902.4545. [43] S. Durr et al. , Science , 1224 (2008), 0906.3599.[44] S. R. Sharpe, Phys. Rev. D46 , 3146 (1992), hep-lat/9205020.[45] C. Amsler et al. (Particle Data Group), Phys. Lett.
B667 , 1 (2008).[46] R. Mawhinney (RBC), PoS
LAT2009 , 081 (2009), 0910.3194.[47] E. E. Scholz (2009), 0911.2191.[48] R. Mawhinney (RBC-UKQCD Collaboration)To be published in the proceedings of the XXVIII International Symposium on Lattice Field Theory.[49] A. Bazavov et al. (The MILC), PoS
LAT2009 , 079 (2009), 0910.3618.[50] C. Aubin et al. , Phys. Rev.
D70 , 094505 (2004), hep-lat/0402030.[51] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa, and A. Vladikas, Nucl. Phys.
B445 , 81 (1995),hep-lat/9411010.[52] Y. Aoki (RBC-UKQCD), PoS
LATTICE2008 , 222 (2008), 0901.2595.[53] Y. Aoki, PoS
LAT2009 , 012 (2009), 1005.2339.[54] J.-R. Cudell, A. Le Yaouanc, and C. Pittori, Phys. Lett.
B454 , 105 (1999), hep-lat/9810058.[55] L. Giusti and A. Vladikas, Phys. Lett.
B488 , 303 (2000), hep-lat/0005026.[56] K. G. Chetyrkin and A. Retey, Nucl. Phys.
B583 , 3 (2000), hep-ph/9910332.[57] J. A. Gracey, Nucl. Phys.
B662 , 247 (2003), hep-ph/0304113.[58] T. van Ritbergen, J. A. M. Vermaseren, and S. A. Larin, Phys. Lett.
B400 , 379 (1997), hep-ph/9701390.[59] M. Constantinou, V. Lubicz, H. Panagopoulos, and F. Stylianou, JHEP , 064 (2009), 0907.0381.[60] P. Di Vecchia and G. Veneziano, Nucl. Phys. B171 , 253 (1980).[61] H. Leutwyler and A. V. Smilga, Phys. Rev.
D46 , 5607 (1992).[62] M. F. Atiyah and S. I. M., Bull. Amer. Math. Soc. , 422 (1963).[63] M. Luscher, Commun. Math. Phys. , 39 (1982).[64] P. de Forcrand, M. Garcia Perez, and I.-O. Stamatescu, Nucl. Phys. B499 , 409 (1997), hep-lat/9701012.[65] M. Falcioni, M. L. Paciello, G. Parisi, and B. Taglienti, Nucl. Phys.
B251 , 624 (1985).[66] M. Albanese et al. (APE), Phys. Lett.
B192 , 163 (1987).[67] Y. Aoki et al. , Phys. Rev.
D73 , 094507 (2006), hep-lat/0508011.[68] Y.-Y. Mao and T.-W. Chiu (TWQCD), Phys. Rev.
D80 , 034502 (2009), 0903.2146.[69] M. Lightman and E. Goode (RBC-UKQCD Collaboration) To be published in the proceedings of the XXVIII International Symposium on Lattice Field Theory.[70] P. Boyle et al. , IBM Journal of Research and Development
49, number 2/3 , 351 (2005).[71] P. A. Boyle, C. Jung, and T. Wettig (QCDOC), ECONF
C0303241 , THIT003 (2003), hep-lat/0306023.[72] P. A. Boyle et al. , J. Phys. Conf. Ser. , 129 (2005).[73] L. Del Debbio, L. Giusti, M. Luscher, R. Petronzio, and N. Tantalo, JHEP , 082 (2007), hep-lat/0701009.[74] J. D. Bratt et al. (LHPC) (2010), 1001.3620.[75] V. Furman and Y. Shamir, Nucl. Phys. B439 , 54 (1995), hep-lat/9405004.[76] N. Christ (RBC and UKQCD), PoS
LAT2005 , 345 (2006).[77] T. Blum et al. , Phys. Rev.
D69 , 074502 (2004), hep-lat/0007038.[78] In this Introduction we combine the statistical and systematic errors in the results. The separate errorsare presented in the following sections.[79] We use the convention that the prime ′ in m ′ res ( m f ) implies that the corresponding residual mass hasbeen determined at a particular value of the light quark mass. m res (without the ′ ) is defined by m res ≡ m ′ res ( ) .[80] Of course the varying quark masses m ud ( b ) and m s ( b ) will also appear in the coefficients of these O ( a ) terms but when expressed in physical units such mass dependence will be of order a .[81] Since N , the number of different b s for which we currently have results, is 2, there is only a single setof ratios R a , Z l and Z h . When specifically discussing our data, we therefore drop the superfix andsimply write R a , Z l and Z h .[82] For completeness, since those have not been included in Appendix C of [1], the finite volume correc-tions for the kaon decay constant in SU ( ) -ChPT read: D L f K xy = − p f " c x + c l d r c x + c l L ! + c l − c x d ( √ c x L ) ..