{'fraction_non_alphanumeric': 0.043198598658931664, 'fraction_numerical': 0.01965456189256385, 'mean_word_length': 4.030557134735646, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Light pseudoscalar decay constants, quark masses, and low energy constants from three-flavor lattice QCD

As part of our program of lattice simulations of three flavor QCD with improved staggered quarks, we have calculated pseudoscalar meson masses and decay constants for a range of valence quark masses and sea quark masses on lattices with lattice spacings of about 0.125 fm and 0.09 fm. We fit the lattice data to forms computed with staggered chiral perturbation theory. Our results provide a sensitive test of the lattice simulations, and especially of the chiral behavior, including the effects of chiral logarithms. We find: f_\pi=129.5(0.9)(3.5)MeV, f_K=156.6(1.0)(3.6)MeV, and f_K/f_\pi=1.210(4)(13), where the errors are statistical and systematic. Following a recent paper by Marciano, our value of f_K/f_\pi implies |V_{us}|=0.2219(26). Further, we obtain m_u/m_d= 0.43(0)(1)(8), where the errors are from statistics, simulation systematics, and electromagnetic effects, respectively. The data can also be used to determine several of the constants of the low energy effective Lagrangian: in particular we find 2L_8-L_5=-0.2(1)(2) 10^{-3} at chiral scale m_\eta. This provides an alternative (though not independent) way of estimating m_u; 2L_8-L_5 is far outside the range that would allow m_u=0. Results for m_s^\msbar, \hat m^\msbar, and m_s/\hat m can be obtained from the same lattice data and chiral fits, and have been presented previously in joint work with the HPQCD and UKQCD collaborations. Using the perturbative mass renormalization reported in that work, we obtain m_u^\msbar=1.7(0)(1)(2)(2)MeV and m_d^\msbar=3.9(0)(1)(4)(2)MeV at scale 2 GeV, with errors from statistics, simulation, perturbation theory, and electromagnetic effects, respectively.


I. INTRODUCTION
Using lattice QCD techniques, the masses and decay constants of light pseudoscalar mesons can be determined with high precision at fixed quark mass and lattice spacing.
Assuming that the chiral and continuum extrapolations are under control, one can therefore calculate from first principles a number of physically important quantities, including • Pion and kaon leptonic decay constants, f π and f K , and their ratio.
• Quark mass ratios, such as m s /m, wherem is the average of the u and d quark masses, and m u /m d .
• Absolute quark mass values, if the mass renormalization constant is known perturbatively or nonperturbatively.
The comparison of f π and f K with experiment provides a sensitive test of lattice methods and algorithms. A precise determination of f K , or f K /f π may in fact be turned around to determine the magnitude of the CKM element V us , as emphasized recently by Marciano [2]. The quark masses are fundamental parameters of the Standard Model, and hence are phenomenologically and intrinsically interesting. Of special importance here is the up quark mass: if m u or m u /m d can be bounded away from zero with small enough errors, it can rule out m u = 0 as a solution to the strong CP problem [3,4]. 1 Finally, the Gasser-Leutwyler parameters give a concise summary of the properties of low energy QCD. In particular the combination 2L 8 − L 5 provides an alternative (although not independent) handle on the up quark mass [5,7].
Extracting these important quantities is predicated on being able to control the chiral and continuum extrapolations. The improved staggered (Kogut-Susskind, KS) quarks [8,9] used here have the advantage in this respect of allowing us to simulate at quite small quark mass: Our lowest m π /m ρ value is ≈ 0.3, a pion mass of roughly 250 MeV. On the other hand, these extrapolations are complicated by the fact that a single staggered quark field describes four species of quarks. We call this degree of freedom "taste" to distinguish it from physical flavor. We simulate the latter by introducing distinct staggered fields for each nondegenerate quark flavor; while we handle the former by taking the fourth root of the staggered quark determinant.
The fact that taste symmetry is violated at finite lattice spacing leads to both practical and theoretical complications. The improvement of the fermion action [8] reduces the splittings among pseudoscalar mesons of various tastes to O(α 2 S a 2 ); yet the splittings are still numerically large, especially on our coarser lattices. This practical problem makes it impossible to fit our data with continuum chiral perturbation theory ( χ PT) expressions (see Refs. [10] and [11], as well as discussion in Sec. IX D). Instead, we must use "staggered chiral perturbation theory" (S χ PT) [12,13,14,15], which includes discretization effects within the chiral expansion. Using S χ PT, we can take the chiral and continuum limits at the same time, and arrive at physical results with rather small systematic errors.
Theoretically, it is not obvious that, in the presence of taste-violations, the fourth root procedure commutes with the limit of lattice spacing a → 0. Assuming that perturbation theory for the standard KS theory without the fourth root correctly reproduces a continuum four-taste theory, then the fourth root trick is correct in perturbation theory [16], since it just multiplies each virtual quark loop by 1/4. However, nonperturbatively, the fourth root version is almost certainly not ultra-local at finite lattice spacing, and the possibility remains that it violates locality (and therefore universality) in the continuum limit. We believe that existing checks [17,18,19,20] of the formalism against experimental results already make this possibility unlikely. The current work adds more evidence that the method gives results that agree well with experiment and have the proper chiral behavior, up to controlled tasteviolating effects that vanish in the continuum limit. However the question is not yet settled.
We discuss this further in Sec. IX D 7 and briefly refer to other recent work that addresses the issue.
This violation of taste symmetry arises because the full axial symmetry (at m q = 0) is broken to a single U(1) subgroup on the lattice. This means that only one of the pseudoscalars, which we call the "Goldstone meson", has its mass and decay constant protected from renormalization. A study of pion masses and decay constants by the JLQCD collaboration [21] explored the masses and decay constants of all of the pseudoscalars in a quenched calculation. We concentrate almost exclusively here on Goldstone mesons, thus avoiding the necessity for renormalization.
We have generated a large "partially quenched" data set of Goldstone meson masses and decay constants using three flavors of improved KS sea quarks. These quantities have been computed with a wide range of sea quark masses (with m u = m d = m s ), and on lattices with lattice spacings of about 0.125 fm and 0.09 fm. We have 8 or 9 different valence quark masses available for each set of sea quark masses and lattice spacing. This data may be fit to chiral-logarithm forms from S χ PT, which at present have been computed for Goldstone mesons only [14,15]. However, since the masses for mesons of other tastes enter into the one-loop chiral logarithms of the Goldstone mesons, some control over those masses is also needed. We have computed most non-Goldstone "full QCD" (valence masses equal to sea masses) pion masses on most of our lattices. We can fit that data to the tree-level (LO) S χ PT form, and use the results for splitting and slopes as input to the NLO terms for the Goldstone mesons. There is, of course, a NNLO error in this procedure, which we estimate in Sec. VI B.
The outline of the rest of this paper is as follows: Section II explains the methodology used to compute raw lattice results (at fixed a and fixed quark mass). In Sec. III, we describe the details of our simulations. We present a first look at the raw data in Sec. IV. Taste violations are discussed in Sec. V, followed in Sec. VI by a detailed description of our S χ PT fitting forms. Relevant results from weak coupling perturbation theory are collected in Sec. VII.
At the current level of precision, electromagnetic and isospin-violating effects cannot be ignored, and we discuss the necessary corrections and the attendant systematic errors in Sec. VIII. Section IX then presents the S χ PT fits, including a description of fit ranges (in quark mass), an inventory of all fit parameters, the resulting fits, and a discussion of various issues relevant to the extraction of physical results. The discussion includes details of the continuum extrapolation, the evidence for chiral logarithms, an estimate of the systematic errors associated with using a (slightly) mass-dependent renormalization scheme, a critical look at the applicability and convergence of the chiral perturbation theory on our data set, bounds on residual finite volume effects, and some comments relevant to the fourth-root trick. In Sec. X, we present our final results, tabulate the systematic errors, and discuss prospects for improving the current determinations.
In collaboration with the HPQCD and UKQCD groups, we have previously reported results for m MS s , the average u-d quark massm MS , and m s /m [22]. The data sets and chiral fits described in detail here are the same ones that were used in Ref. [22].

II. METHODOLOGY
For the axial current corresponding to the unbroken (except by quark mass) axial symmetry, the decay constant f P S can be found from the matrix element ofψγ 5 ψ between the vacuum and the pseudoscalar meson. In terms of the one component staggered fermion field ψγ 5 ψ corresponds to the operator O P (t) =χ a ( x, t)(−1) x+t χ a ( x, t) . (1) Here a is a summed color index. The relevant matrix element can be obtained from a pseudoscalar propagator using O P as both the source and sink operator: where m P S is the mass of the pseudoscalar, and V s is the spatial volume.
The decay constant is obtained from C P P by [21,23] f P S = (m x + m y ) V s 4 where m x and m y are the two valence quark masses in the pseudoscalar meson. Throughout this paper we use the convention where the experimental value of f π is approximately 131 MeV. Note that in computing this meson propagator we must take care to normalize the lattice Dirac matrix as M = am + D /. The four in the denominator arises from the number of tastes natural to the Kogut-Susskind formulation. (See unnumbered equations between Eqs. 7.2 and 7.3 in Ref. [23].) 2 However, the point operator O P has large overlap with excited states. For calculating masses it is customary to use an extended source operator that suppresses these overlaps, together with a point sink. In our case, this extended operator is a "Coulomb wall," i.e., we fix to the lattice Coulomb gauge and sum over all lattice points on a timeslice: x, yχ ( x, t)(−1) x+t χ( y, t) .
We can calculate propagators with any source or sink operator we wish. Ignoring excited state contributions, we have for example We will use the shorthand "PP" for point-source point-sink propagators, "WP" for Coulombwall-source point-sink propagators, "PW" for point-source Coulomb-wall-sink propagators, and "WW" for Coulomb-wall source and sink propagators. In previous calculations of pseudoscalar decay constants the relation C P P = C 2 W P /C W W has often been used to get the point-point amplitude. However, the wall-wall propagator has large statistical fluctuations and severe problems with excited states, as was discussed in Ref. [21]. To be able to use the PP operator to get C P P directly, rather than indirectly by way of the ratio formula, one needs much better statistics. We do this by replacing the point source with a "random-wall" source, which simulates many point sources. We set the source on each site of a time slice to a three component complex unit vector with a random direction in color space, and use this as the source for a conjugate gradient inversion to compute the quark propagator, whose magnitude is squared to produce the Goldstone pion propagator. Thus, contributions to a meson propagator where the quark and antiquark originate on different spatial sites will average to zero and, after dividing by the spatial lattice volume, this source can be used instead of O P .

Figures 1 and 2 show masses and amplitudes from pion propagators with random-wall
and Coulomb-wall sources and sinks. In Fig. 1, we can see that extraction of masses from the "WW" propagators is almost hopeless. Including an excited state helps, but statistical errors become very large. In Fig. 2, the WW amplitudes are also slower to plateau, though not as bad as the masses. As a consistency check, note that the WP and PW amplitudes are equal, and the masses extracted from the diagonal PP and WW propagators approach their value from above (since excited states must contribute to these propagators with the same sign as the ground state). As an additional illustration of the difficulties with using the Coulombwall-Coulomb-wall propagator, Fig. 3 plots the ratio of the point-point pion propagator (using the random wall source) to the alternative P P W P W P /P W W (with a different mass than in Figs. 1 and 2). While this ratio is approaching one, it is clear that we would either need very large minimum time in the fit or a careful removal of excited states to use the "WW" propagators.
Given the problems with the WW propagators, we have opted to use only the Coulombwall-point-sink and random-wall-point-sink propagators. We performed a simultaneous fit to these two propagators, with an amplitude for each propagator and a common mass. In these fits the WP propagator dominates the determination of the mass; while the amplitude of the PP propagator is required for computing the decay constant. Since the combination C P P /m 3 π is needed for determining f π and the mass and amplitude in a fit to a meson propagator are strongly correlated, we used this combination as one of our fitting parameters.
That is, we fit the point-point and wall-point meson correlators to so that A P P is the desired combination C P P /m 3 π . Since the correlation between m π and the propagator amplitude is positive, the statistical error on the quantity C P P /m 3 π is somewhat smaller than a naive combination of the errors on C P P and m π .

III. SIMULATIONS
These calculations were made on lattices generated with a one loop Symanzik and tadpole improved gauge action [9,24] and an order a 2 tadpole improved Kogut-Susskind quark action [8]. Parameters of most of the lattices, as well as the light hadron spectrum, are in Ref. [10,19]. The determination of the static quark potential, used here to set the lattice spacing, is presented in Ref. [10,19,25]. In addition to the runs tabulated in Ref. [  The masses shown here come from a separate spectrum calculation, using more source time slices than were used in the partially quenched calculations, and using more lattices at am ′ = 0.03.

Equation 7
can be used to express these masses and decay constants in units of r 1 . about 2.5 fm with the exception of the am ′ , am ′ s = 0.005, 0.05 run, where the spatial size is about 3.0 fm. Table I lists the parameters of the runs used here.
We note here that the values for am ′ s were approximately tuned from the vector to pseudoscalar meson mass ratio in initial runs with fairly heavy quarks. Our best determinations of the physical strange quark mass at these lattice spacings turned out to be lower by 8 to 22% (coarse) and 6 to 12% (fine) than the nominal values m ′ s , where the range depends on whether or not taste-violating terms (as determined by S χ PT fits) are set to zero before demanding that m π and m K take their physical values on a given lattice. Pseudoscalar propagators were calculated on lattices separated by six units of simulation time, using two source time slices per lattice. For the coarse lattices, nine valence quark masses were used, ranging from 0.1m ′ s to m ′ s ; while for the fine lattices eight masses ranging from 0.14m ′ s to m ′ s were used. In all but one of the runs, the source slices were taken at different points in successive lattices, which leads to smaller autocorrelations than using the same source time slices on all lattices. The effects of the remaining correlations among the sample lattices were estimated in two ways. First, jackknife error estimates for the masses am ′ / am ′ s 10/g 2 ∆m 2 (4) (4) ∆f (1)   Fitting the pseudoscalar propagators produces masses and decay constants in units of the lattice spacing a, and to convert to physical units we must estimate a from a calculation of some dimensional quantity whose value is known. This amounts to saying that we are calculating ratios of these quantities to some other quantity calculated from these simulations.
We express our results in units of a length obtained from the static quark potential, r 1 , where r 2 1 F (r 1 ) = 1.0 [25,26]. This has the advantage that r 1 can be accurately determined in units of the lattice spacing. But r 1 is not a directly measurable quantity, and its physical value must in turn be obtained from some other quantities. We have calculated the static quark potential in all of these runs, and fit it to determine r 1 /a. To smooth out statistical fluctuations in these values, we then computed a "smoothed r 1 " by fitting the r 1 /a values to a smooth function. A simple form, which gives a good fit over the range of quark masses and gauge coupling used here, is [19] log(r 1 /a) = C 00 + C 10 (10/g 2 − 7) + C 01 am tot + C 20 (10/g 2 − 7) 2 , where m tot = 2m ′ + m ′ s . The results of the fit are C 00 = 1.258(3) C 10 = 0.937(9) When we need an absolute lattice scale, we start with the scale from Υ 2S-1S or 1P -1S splittings, determined by the HPQCD group [17,27]. This gives a scale a −1 = 1.588 (19) GeV on the coarse 0.01/0.05 lattices, and a −1 = 2.271(28) GeV on the fine 0.0062/0.031 lattices. For light quark masses < ∼ m s /2, the mass dependence of these quantities and of r 1 appears to be slight, and we neglect it. With our smoothed values of r 1 /a, we then get r 1 = 0.324(4) fm on the coarse lattices and r 1 = 0.320(4) fm on the fine lattices.
To extrapolate r 1 to the continuum, we first assume that the dominant discretization errors go like α S a 2 . Using α V (q * ) [28] (with scale q * = π/a) for α S gives a ratio (α S a 2 ) fine /(α S a 2 ) coarse = 0.427. Extrapolating away the discretization errors linearly then results in r 1 = 0.317(7) fm in the continuum. However, taste-violating effects, while formally O(α 2 S a 2 ) and hence subleading, are known to be at least as important as the leading errors in some case. Therefore, one should check if the result changes when the errors are assumed to go like α 2 S a 2 . Taking α S = α V (3.33/a) gives a ratio (α 2 S a 2 ) fine /(α 2 S a 2 ) coarse = 0.375; while a direct lattice measurement of the taste-splittings gives a ratio of 0.35. Extrapolating linearly to the continuum then implies r 1 = 0.318(7) fm or r 1 = 0.319(6) fm, respectively, in agreement with the previous result. For our final result, we use an "average" ratio of 0.4 and add the effect of varying this ratio in quadrature with the statistical error. We obtain r 1 = 0.317(7) fm. A systematic error of 0.03 fm in r 1 from our choice of fitting methods is omitted since it is common to all our runs and cancels out in the final results here. Using our current value r 0 /r 1 = 1.472(7), the result for r 1 implies r 0 is about 7% smaller than the standard phenomenological choice r 0 = 0.5 fm, although the difference is within the expected range of error of the phenomenological estimates [26]. These figures also show a crude extrapolation to a = 0, made by taking a linear extrapolation in α S a 2 using pairs of points with the samem ′ /m ′ s . In Fig. 7 one pair of extrapolated points has diagonal lines showing the data points that were extrapolated to produce this point. In hindsight, m ′ s used in the a ≈ 0.09 fm runs was smaller than that used in the a ≈ 0.125 fm runs, as indicated by the fact that the finer lattice points fall slightly to the left of the corresponding coarse lattice points.

V. TASTE SYMMETRY VIOLATIONS
As mentioned above, we use the term "taste" to denote the different staggered-fermion species resulting from doubling. At finite lattice spacing, taste symmetry is violated. Although the improved staggered action reduces the taste violating effects to O(α 2 S a 2 ) from O(α S a 2 ) with unimproved staggered fermions, the violations are still quite significant numerically. Figure 9 shows the splittings between pions of various tastes on our coarse lattices. There are 16 such pions, π B , where B = 5, µ5, µν, µ, I (ν > µ) labels taste matrices in the taste Clifford algebra generated by Euclidean gamma matrices ξ µ . The π 5 is the Goldstone (pseudoscalar taste) pion, whose mass is required to vanish in the chiral limit by the exact (nonsinglet) lattice axial symmetry. All the pions in Fig. 9 are flavor-charged, i.e., π + mesons.
Thus there are no contributions from disconnected graphs, even for the taste-singlet π + I . The approximate "accidental" SO(4) identified by Lee and Sharpe [12] is clearly a good symmetry: there is near degeneracy between π + 05 and π + i5 , between π + 0i and π + ij , and between π + 0 and π + i . When we assume such degeneracy, we can think of the index B as running over the multiplets 5, A, T, V, I with degeneracies 1, 4, 6, 4, 1, respectively.
The fit in Fig. 9 is to the tree-level chiral form given in Refs. [12,14]: The slope, µ tree , is the same for all tastes, but there are constant splittings for each non-  Goldstone multiplet (∆ 5 = 0). Although the fit is poor (chiral logs, including tasteviolations, are needed), it does give the pion squared masses within a few per cent: The biggest deviation, 7%, is for the Goldstone pion at the lowest mass; most other deviations are ∼ 2%. Table III shows the values of a 2 ∆ B coming from the fit on the coarse lattices. On the fine lattices, we have measured non-Goldstone pion masses only on the set with quark masses 0.0124, 0.031. So we directly compare the splittings with those of the corresponding coarse lattice (masses 0.02, 0.05). The fine-lattice splittings are smaller by a common factor of 0.35, within errors. This is consistent with the expectation that taste violations go like O(α 2 S a 2 ). Indeed, if we take α S = α V (q * ) [28] and choose q * = π/a because taste violations occur at the scale of the cutoff, we find The ratio of taste-violating terms between fine and coarse lattices is an input to the chiral fits for Goldstone pions discussed below. The measured splitting ratio of 0.35 is used as a central value. The error can be estimated by varying q * in Eq. (10): q * = π/(2a) gives a ratio of 0.324; while 2π/a gives 0.398. We take 0.3-0.4 as an appropriate range for our analysis of systematics.
We warn the reader here that the notation in Eq. (9) can be slightly misleading. We have shown explicitly the a 2 factor in the taste-violating splitting, a 2 ∆ B , but this does not mean that ∆ B itself is independent of lattice spacing, or even that it approaches a non-zero constant in the continuum limit. Indeed, the argument above implies that ∆ B is a slowly varying function of a that goes like α 2 V (π/a) for small a. A similar comment applies to the other taste-violating parameters introduced in Sec. VI A: the a 2 dependence is always shown explicitly, but dependence on a through the coupling is hidden.
In physical units, the splittings on the coarse lattices are quite large. The largest is for the taste-singlet pion: a 2 ∆ I ≈ (450 MeV) 2 ; while the smallest, for the taste axial-vector pion, is a 2 ∆ A ≈ (280 MeV) 2 . Given the size of these splittings, which are discretization errors, it is not surprising that the lattice data is not well fit by continuum chiral perturbation theory ( χ PT) forms. Figure 10 shows such an attempted fit for the Goldstone f π to the standard NLO partially quenched continuum form [29] plus analytic NNLO terms. More details about this fit will be explained below, when we discuss the corresponding fits that take into account taste violations. For the moment, we simply remark that the minuscule confidence level (CL ≈ 10 −250 ; χ 2 /d.o.f. = 8.77 with 204 degrees of freedom) shows how hard it is to ignore lattice artifacts at the level of chiral logarithms.

VI. STAGGERED CHIRAL PERTURBATION THEORY
Lee and Sharpe [12] found the chiral Lagrangian that describes a single staggered field.
Their Lagrangian includes the effects of taste violations at O(a 2 ) as well as the standard violations of chiral symmetry from mass terms at O(m q ), where m q is a generic quark mass.
They introduced a power counting that considers m q and a 2 to be of the same order, which is appropriate here: In Fig. 9 the splittings are comparable to the squared meson masses. Tree level (LO) is thus O(m q , a 2 ); chiral logs appear at one-loop (NLO) and are O(m 2 q , m q a 2 , a 4 ). 3 The Lee-Sharpe Lagrangian is not directly appropriate to the calculations here because it has only one flavor (one staggered field). Aubin and Bernard [13,14,15] have generalized Ref. [12] to n staggered flavors and shown how to accommodate the 4 √ Det trick in loop calculations. This is what is meant by "staggered chiral perturbation theory," S χ PT.
Continuum chiral perturbation theory can be thought of as an expansion in the dimen- 3 Throughout this paper, we define the order of a contribution to be the order of the corresponding term in the chiral Lagrangian. This is the simplest way to keep the power counting consistent between decay constants and meson masses, although it does lead to the unnatural statement that the tree-level f π is "O(m q )" since it comes from the kinetic energy term in the chiral Lagrangian. What matters ultimately is only the relative size of contributions: the first correction to the tree-level value of m 2 π or f π is smaller by one power of m q . sionless quantity where 2µm q is the tree-level mass of a qq meson. For physical kaons, we expect the relevant quark mass parameter to be χ ud,s ≡ (χ ud + χ s )/2 ≈ 0.18 (where χ ud is the average value for the u and d quarks); this is reasonable given the experimental result f K /f π ≈ 1.22.
Staggered chiral perturbation theory is a joint expansion in χ q and χ a 2 , which measures the size of the O(a 2 ) taste violations: where a 2 ∆ is a "typical" taste-violating term. Taking for a 2 ∆ the average meson splitting  [14,15], a 2 δ ′ A , to estimate ∆ and χ a 2 , one gets slightly smaller values.

A. NLO forms
One-loop chiral logs and analytic terms have been calculated in S χ PT for Goldstone meson masses [14] and decay constants [15]. Partially quenched results are included, so all forms needed to fit the numerical data are available.
References [14,15] express the chiral logarithms in terms of where Λ χ is the chiral scale, and L is the spatial dimension. The finite volume correction terms δ 1 and δ 3 are [13] δ 1 (mL) = 4 r =0 where K 0 and K 1 are Bessel functions of imaginary argument, and r, which labels the various periodic images, is a three-dimensional vector with integer components. We have assumed here that corrections due to the finite time extent are negligible; this is true for our lattices, for which the time dimension is between 2.7 and 3.4 times greater than the spatial dimension.
The function ℓ(m 2 ) in Eq. (13) arises from tadpole diagrams with a single meson propagator; ℓ(m 2 ) in Eq. (14) comes from double-poles, which are present only in the partially quenched (and quenched) cases, not in the full QCD limit. In practice, we compute the sum in Eq. (15) or Eq. (16) with cutoff | r| ≤ N, where N is an integer, and increment N by 1 until the sum changes by a fractional amount ≤ ǫ. To be conservative, we take ǫ = 10 −9 for central-value fits. However, a much weaker criterion, ǫ = 0.001, is adequate to reduce the error in the sum well below our statistical errors, and we often use the weaker criterion for alternative fits in the systematic error estimates.
In the generic case relevant to our data (m u = m d ≡m = m s and no degeneracies between valence and sea quarks), the NLO S χ PT expressions for a meson P composed of valence quarks x and y are [14,15] (m NLO Here µ and f are the continuum chiral parameters, δ ′ V and δ ′ A are LO taste-violating parameters (hairpins), L i are the NLO Gasser-Leutwyler [1] coefficients, and L ′ and L ′′ are linear combinations of the taste-violating NLO coefficients. The reason for using the tree-level µ tree parameter from Eq. (9) in the L i terms will be explained in Sec. VI B. X Ξ and Y Ξ are flavor-neutral mesons of taste Ξ made of x,x and y,ȳ quarks, respectively, and U Ξ , D Ξ , and S Ξ are corresponding flavor-neutral mesons made from u, d, and s sea quarks, respectively.
The index Q runs over the 6 mesons made from one valence and one sea quark, and B runs over the 16 meson tastes. The residues R Refs. [14,15]. For completeness, we quote them here: .
Each of these residues is a function of two sets of masses, the "denominator" set {M} = {m 1 , m 2 , . . . , m n } and the "numerator" set {µ} = {µ 1 , µ 2 , . . . , µ k }. The indices j and i, 1 ≤ j, i ≤ n, refer to particular denominator masses; the prime on the product in the denominator of Eq. (19) means that ℓ = j is omitted.
In Eqs. (17) and (18), the denominator mass-set arguments are shown explicitly; they are The index j in Eqs. (17) and (18) is summed over the denominator masses. Sets for axialvector taste (A) are found from the corresponding vector taste (V ) sets by taking V → A in Eq. (21). The masses m η I , m η V , m η ′ V are given by [14] .
The numerator mass-set arguments of the residues in Eqs. (17) and (18) are not shown explicitly because they are always where the taste label Ξ is taken equal to the taste of the denominator set.
Degeneracies among the various masses in Eqs. (17) and (18)  . Formulas for many of these degenerate cases appear in Refs. [14,15]. For the other cases, one can carefully take appropriate limits in Eqs. (17) and (18), or, more conveniently, return to the original integrands in Refs. [14,15] and take the limits before performing the momentum integrations.
Because Eqs. (17) and (18) are quite complicated, it is useful to write down a simple result that shows more clearly how taste violations change the continuum chiral behavior.
The pion decay constant in full QCD with two (degenerate) flavors (m x = m y = m u = m d ≡ m, with the strange quark integrated out) is particularly simple. In that case, the result corresponding to Eq. (18) is with B running as usual over the 16 possible tastes, and In Eq. (24), the term multiplied by −2 gives the average of all tastes and becomes the standard SU(2) L × SU(2) R chiral logarithm in the continuum limit, when all tastes are degenerate. The terms multiplied by −4 clearly vanish in the continuum limit because Sincem ′ is significantly less than m ′ s for many of our runs, Eq. (24) is often not a bad approximation to the chiral behavior of our (full QCD) data. It will be useful in the discussion of finite volume effects in Sec. IX D 6.
We note here that Refs. [12,14,15] explicitly include in the chiral Lagrangian the effects of terms in the O(a 2 ) staggered-quark Symanzik action that violate the taste symmetries.
There are also "generic" O(a 2 ) terms in the Symanzik action that have the same symmetries as the continuum QCD action and are not included explicitly. An example is a 2ψ D 2 D µ (γ µ ⊗ I)ψ, where γ µ and I act on spin and taste indices, respectively. The effect of such terms on the chiral Lagrangian is to produce O(a 2 ) variation in physical parameters such as f , µ, and, at higher order in m q , L i . We build the possibility of such generic variation in physical parameters into the chiral fits below. 4 Since our staggered action is a 2 tadpole improved, we expect such generic variation to be of size α S a 2 Λ QCD ≈ 2%. When we extrapolate the physical parameters to the continuum, we will need to know how α S a 2 changes from the coarse to fine lattice. As in the case of taste violations, such discretization errors occur at the scale of the cutoff. Therefore, we use α S = α V (q * = π/a) for central values, and allow q * to vary between π/(2a) and 2π/a for the error estimate. We have and a range for this ratio of 0.398 to 0.441.
As they stand, Eqs. (17) and (18) are slightly inconvenient because the renormalization of the O(m q a 2 ) analytic NLO parameters L ′ and L ′′ under a change in the chiral scale Λ χ is complicated and involves the physical L i parameters. This is due to the fact that the meson masses multiplying the logarithms include O(a 2 ) splittings. It is more natural, therefore, to redefine the L ′ and L ′′ by associating particular O(a 2 ) terms with the L i . We make the replacements 16µ tree where a 2 ∆ I is given by Eq. (9), and is the average splitting. After these redefinitions, a change in Λ χ renormalizes L ′ according to: while L ′′ is independent of scale. From this we would expect that L ′ is comparable in size to δ ′ A + δ ′ V , an expectation that is borne out by the fits. The L i renormalize by: with

B. NNLO Terms
As we will see below, the high statistical precision of our data requires us to go beyond the NLO formulas, even for subsets of the data that include only the lighter valence quark masses. We include explicitly all NNLO physical analytic parameters, i.e., all analytic terms of O(m 3 q ). There are five such terms for m 2 P + and an additional five for f P + [11]. Expressed in terms of χ q defined in Eq. (11), they are given by: where "NLO" denotes the lower order contributions, (the corrections to the leading "1" in Eqs. (17) and (18), with the substitutions in Eq. (27)). The interchange symmetries among valence quarks x ↔ y and sea quarks u ↔ d ↔ s ↔ u restrict the form of the NNLO corrections. These terms were obtained independently in Ref. [30].
Possible analytic taste-violating terms at NNLO of O(m 2 q a 2 ) are included implicitly by allowing the L i (O(m 2 q ) terms) to vary with lattice spacing. 5 However, such variation can be caused either by taste-violating terms in the Symanzik Lagrangian or by terms with the same symmetries as the continuum operators but with explicit factors of a 2 (i.e., by generic discretization errors on the L i ). As explained above, the generic discretization errors are expected to be ≈ 2%; while the new taste-violating terms could change the apparent value of L i between coarse and fine sets by order χ coarse a 2 − χ fine a 2 ≈ 6%. For our preferred fits, we use Bayesian priors [31] to restrict the differences in the L i on coarse and fine sets to be at most 7.5% (for a 3-σ variation); while in alternative fits used to assess the systematic error, we relax or tighten this restriction. (See Sec. IX B.) Note that when we extrapolate the L i to the continuum, we have no a priori way to distinguish variation as α S a 2 (generic discretization errors) from variation as α 2 S a 2 (taste violations). Therefore, we consider both types of variation (i.e., Eqs. (26) and (10), as well as the ranges of these ratios) and include the difference in the systematic error. In practice, these alternative fits and these assumptions in how the L i are extrapolated to the continuum contribute only a small fraction of the total systematic error. (See Table VII.) At NNLO, O(m q a 4 ) analytic terms may also affect m 2 P + 5 and f P + 5 . We neglect such terms, which would make contributions similar to that of L ′ and L ′′ in Eqs. (17) and (18), but multiplied by a 4 instead of a 2 . Because χ a 2 is only 0.09 in the worst case, we would expect the error thereby induced to be at most χ 2 a 2 ∼ 1%. In fact, since we consider a generous range for how taste-violating a 2 terms may vary as we go from coarse to fine lattices (see discussion immediately following Eq. (10)), the effects of O(m q a 4 ) analytic terms should already be included in our systematic error estimates.
The final possible NNLO terms in the Lagrangian are O(a 6 ). However, it is easy to see that such terms do not contribute to m 2 P + 5 and f P + 5 . The Goldstone theorem requires that m 2 P + 5 be proportional to at least a single power of quark mass (m x + m y in this case), and terms in the Lagrangian must have at least two derivatives to make analytic contributions to f P + 5 through Noether's theorem or wave-function renormalization. In addition to analytic terms at NNLO, there are, of course, NNLO chiral logarithms (from 2-loop graphs, as well as 1-loop graphs that involve NLO parameters). These nonanalytic terms have not been calculated in S χ PT, but in any case are not expected to be important here: Wherever the quark masses or splittings are large enough for the analytic NNLO terms to be significant, the NNLO logarithms should be slowly varying and well approximated by analytic terms. As discussed in Sec. IX B, the NNLO terms make a difference primarily in the interpolation around m s , not in the extrapolation tom. The systematic errors inherent in our treatment of the NNLO terms are estimated by varying the masses we fit to and the Bayesian priors governing these terms and their changes with a, as well as by adding still higher (NNNLO) terms.
There are also NNLO effects induced by the ambiguity in the parameters one puts into NLO expressions. In particular, we have at present expressed the "chiral coupling," 1/(16π 2 f 2 ) in Eqs. (17) and (18), in terms of bare (tree-level) parameter f . Replacing f with the experimental value of f π , say, would generate a difference at NNLO. As we discuss below, the difference between f π and f is significant: ≈ 13%. If we had the full NNLO expression, including 2-loop effects, then the ambiguity would be resolved up to terms of NNNLO. But in the present case there is no a priori way to decide this issue.
We argue, however, that putting a physical parameter in the chiral coupling (f → f π ) is likely to result in a better convergent χ PT. This is similar to the argument for using a physical, rather than bare, coupling in weak-coupling perturbation theory [24]. In practice, we consider three versions of the fits: (1) Fix coupling as 1/(16π 2 f 2 π ).
(3) Write coupling as ω/(16π 2 f 2 π ) and treat ω as an additional fit parameter: either allow it to vary freely or force it to vary around 1 using Bayesian priors.
Good fits are possible with all three choices. Both because of the argument above, and because it guarantees that the NLO chiral logarithms for very light quarks have the expected weight relative to the tree level terms, we take choice (1) for our central values. Choice (3), with its extra parameter, results in the highest confidence levels of the three. When ω is allowed to vary freely, its value decreases as higher quark masses are included in the fits, reaching ω ≈ 0.6 by set II (see Sec. IX A). This is similar to replacing f → f K , perhaps not surprising for fits that must cover a range of valence masses up to a large fraction of m s . But for the quantities computed here, all of which are sensitive to the chiral behavior at low quark mass, we do not include fits with ω free since we expect 1/(16π 2 f 2 π ), not 1/(16π 2 f 2 K ), to be the correct weight for the logarithms in the low mass regime. We still allow fits with range ω = 1.0 ± 0.1 because 10% is roughly the difference between the physical f 2 π and its value in the chiral limit. As discussed in Sec. X, the main effects of including fits with arbitrary ω would be to increase the systematic error in f π by about 1 MeV (with a corresponding effect on f K /f π ) and to double the simulation error on m u /m d (which error, however, is small compared to unknown electromagnetic effects).
Good fits with choice (2) require f > ∼ f K (equivalent to ω ≈ 0.6 in choice (3)) and quite large NNLO terms. In addition, the O(a 2 ) NLO parameter a 2 L ′′ becomes unreasonably large (∼ (630 MeV) 2 on the coarse lattices). For these reasons we exclude choice (2) fits from the systematic error estimates; including them would increase systematic errors by amounts comparable to those of arbitrary ω fits.
Similar considerations apply to the parameter µ. In the analytic terms involving the L i in Eqs. (17) and (18) (and hence in Eqs. (32) and (33)), we argue that it is best to put in µ tree from the linear (tree-level) fits, Eq. (9). The L i are then multiplied by actual squared meson masses (within the errors of Eq. (9)). This corresponds to how such terms are interpreted in continuum χ PT analysis (see, e.g., Ref. [1]). An alternative, a posteriori, choice would be to use the chiral limit of m 2 π /(2m) coming from the full NNLO fit. This would replace µ tree in Eqs. (17) and (18) by a number 5% smaller. Since Eq. (9) gives a maximum of 7% errors, we choose that larger value as the systematic effect. It would however be unreasonable to replace µ tree in Eqs. (17) and (18) or (Eqs. (32) and (33)) by the fit parameter µ itself. That is because the effective value of µ is corrected by terms involving m s that do not go away in the chiral limit for the light quarks. Indeed, one might expect corrections of at least In practice the fit parameter µ from our preferred NNLO fit on the intermediate valence mass set (subset II , Sec. IX A) is 29% less than µ tree , which means that µ(m x + m y ) is significantly less than our measured value of m 2 This difference improves to 13% in the fit to the lightest masses.
As discussed in the introduction, all meson masses appearing in the NLO chiral logarithms in Eqs. (17) and (18) are similarly evaluated using the previously determined values of the taste splittings, a 2 ∆ B , and µ tree from the fit of our "full QCD" data for all meson tastes to Eq. (9). In our results for masses and decay constants, the NNLO error introduced by this procedure is negligible. That is because the effect of the small errors in masses in the chiral logarithms on our extrapolated values is almost completely canceled by the effect of variations of the analytic parameters in the fit. We can check this by replacing µ tree in the fit by the (5% different) chiral limit value; the effect is about 0.2% on quark masses and less than 0.1% on the decay constants, in both cases much less than the total systematic error. For the L i , changing the value of µ in the chiral logarithms does not completely cancel the effect of changing its value in the analytic terms, but there is some cancellation.
Therefore the 7% systematic effect in the L i discussed in the previous paragraph remains a conservative estimate of the error.
Since there are, of course, many additional NNNLO terms, it is nonsystematic to include only one each for m P + 5 and f P + 5 . However, we pick these terms involving valence masses because there is a steeper dependence on the valence masses than on the sea quark masses.
For lower quark masses, where we expect χ PT to work well, we fit to Eqs. (34) and (35) only to estimate systematic errors due to the truncation of χ PT. When the fits include valence masses equal to or greater than m ′ s , we also use Eqs. (34) and (35) in order to improve the interpolation around the strange quark mass. In the former case, we find that the values of |ρ (f ) | and |ρ (m) | coming from the fits are typically less than 0.1; in both cases they are always less than 0.2 (including when we fit to Eqs. (36) and (37) -see Table IV).
Another form, used only for interpolations around the strange quark mass, adds on the square of the NLO term as a mock-up of the effect of 2-loop chiral logarithms: where again NLO and NNLO represent the contributions from Eqs. (17), (18), (32) and (33). The absolute values of the new coefficients σ (m) and σ (f ) in the fits are never greater than 0.14.

VII. PERTURBATION THEORY
Because the axial current we use to compute decay constants is partially conserved, there is no need for current renormalization. Mass renormalization is however needed to find continuum (MS) quark masses, as discussed in Ref. [22]. Let Z m be the mass renormalization factor that connects the bare lattice mass (am) 0 and the MS mass at scale Λ: 6 Here, unlike in Ref. [22], we have shown explicitly the plaquette tadpole improvement factor u 0P , necessary because the MILC improved staggered action defines the lattice quark mass in a somewhat unconventional manner.
The renormalization factor Z m enters the analysis in another way. As mentioned above, we need to renormalize the parameter µ if we wish to compare its values at different lattice spacings. More precisely, we need the ratio R m is in principle independent of Λ, although when Z m is evaluated at any given order in perturbation theory, there is residual Λ dependence from neglected higher order terms. For definiteness, we take Λ = 2 GeV. Z m is given by [22] Z m (aΛ) where α V is determined from small Wilson loops using third order perturbation theory [28,32], the optimal scale q * is estimated using a second order BLM method [33], and Z (2) m is [34,35,36] 6 Λ was called µ in [22], a notation we avoid here for obvious reasons.
with the above value of R m defines what we mean by "equality" of the parameter µ on coarse and fine lattices. Of course, Eq. (42) may be violated by generic O(a 2 ) scaling violations (∼ 2%), as well as by perturbative errors. A priori, one expects a two-loop correction to R m of order α 2 V . This is ≈ 6%. In practice, fits have a confidence level that is higher than those of our preferred fits if we take R m ≈ 0.87 to 0.89, i.e., a 7 to 9% difference from R m = 0.958. Although it is not possible to separate the perturbative errors from the discretization errors in this difference, here and in Ref. [22] we take the larger value, 9%, as the conservative estimate of perturbative errors. This is ≈ 1.5α 2 V . For quantities that do not directly involve perturbation theory, such as the decay constants and the ratio m u /m d , we do not quote perturbative errors, per se. But R m still enters the chiral fits, so we include fits with R m = 0.87 among the alternatives.
Another rough estimate of R m comes from µ tree , Eq. (9). Without the proliferation of parameters at NLO and NNLO, the tree-level form makes possible well-controlled fits on coarse and fine lattices separately. We get R m ≈ 0.977. But note that Eq. (9) can have up to ∼ 7% errors in describing the data, and there are also discretization errors in this estimate.

VIII. ELECTROMAGNETIC AND ISOSPIN-VIOLATING EFFECTS
Given the precision we are aiming at here, it is necessary to take into account electromagnetic (EM) and isospin-violating effects, at least in an approximate way. Our simulation is in isospin-symmetric QCD, with m u set equal to m d , and the electromagnetic coupling, e 2 = 4πα EM , set to 0. This means that when we compare meson masses to experiment to determine the physical quark massesm and m s , we must first adjust the experimental numbers to what they would be in a world without EM effects or isospin violation. This is particularly important for the pion, since the difference between m 2 π + and m 2 π 0 is almost 7%. Because the adjustment is only approximate, there are some residual systematic errors on the quark masses, as discussed in Ref. [22].
The decay constants, as well as the low energy constants L i , are by definition pure QCD quantities, so we do not have to take EM effects directly into account in our determination. 7 Nevertheless, there are indirect EM effects on f π and f K , which come in through the quark 7 However, the EM corrections must be explicitly evaluated when the decay constants are compared to experiment [37,38]. masses when we extrapolate to the physical point. Isospin violations are irrelevant for the L i , which are defined to be mass independent. But for the decay constants, there are both direct and indirect isospin-violating effects, which we estimate below. The end result is that both the (indirect) EM and isospin-violating errors on decay constants are very small, as long as we are careful to extrapolate to the appropriate values of the quark mass in each case. However the EM error on m u /m d is large unless we are willing to assume that the EM effects on meson masses are accurately known.
Electromagnetism can be included in χ PT in a systematic way. Dashen's theorem [39] summarizes the EM effects on meson masses at lowest nontrivial order in e 2 and the quark masses. It states that m 2 π + and m 2 K + receive equal O(e 2 ) contributions in the chiral limit; while the π 0 and K 0 masses are unaffected. However, there can be large and different EM contributions to m 2 π + and m 2 K + of order e 2 χ ud,s [40,41,42,43]. Following Ref. [44], we let ∆ E parameterize violations of Dashen's theorem: Then Refs. [40,41,42,43] suggest ∆ E ≈ 1. Most of these corrections are probably to the charged meson masses. Indeed, the violation of Dashen's theorem for the π 0 is O(e 2 χ ud ) [41] and therefore small. The EM contribution to m 2 K 0 , on the other hand, is in principle the same order as the violations of Dashen's theorem for the charged masses, e 2 χ ud,s [41].
Nevertheless, a large N c , extended NJL model calculation [42] finds a tiny EM correction to the K 0 mass at this order. To be conservative, though, we allow for EM contributions to m 2 K 0 of order of half the violations of Dashen's theorem, with unknown sign: The effects of isospin violation in the pion masses are quite small. When m u = m d , m 2 π 0 gets a contribution of order (χ u − χ d ) 2 . The isospin-violating splitting (m π + − m π 0 ) QCD is estimated as 0.17(3) MeV in Ref. [1], and as 0.32 (20) MeV in Ref. [45]. In the kaon system, on the other hand, the effects of isospin violation are clearly important, as is obvious from the fact that the experimental K + -K 0 splitting is of opposite sign to that in Eq. (43) for any ∆ E > −1. In our calculation, we can reduce the isospin violating effects in the kaon masses to the same order as in the pion system by focusing on the isospin averaged quantity (m 2 K 0 + m 2 K + )/2. We then neglect the remaining isospin violations in the meson masses. We have checked, using the estimates for (m π + − m π 0 ) QCD above, that the indirect effect of such isospin violations on decay constants is extremely small: < ∼ 0.03%. These isospin violations were also neglected in the computation of quark masses in Ref. [22]. We note, however, that including isospin violations could have some small effect there, in particular on the result for the ratio m s /m. If (m π + − m π 0 ) QCD is at the upper end of the range in [45], the central value for m s /m in Ref. [22] could be changed from 27.4 to as low as 27.2.
Based on the above discussion, we may determine the physical values ofm and m s by extrapolating the lattice squared meson masses to m 2 π and m 2K , given by where experimental values are to be used on the right hand side. Allowing for EM corrections to the K 0 mass, Eq. (44), replaces ∆ E ≈ 1 in Eq. (45) with an effective value in the range 0-2, which is in any case a conservative range for ∆ E that includes the Dashen theorem result. Here and in Ref. [22], we take ∆ E = 1 for the central value, and use 0 ≤ ∆ E ≤ 2 to estimate systematic errors inm, m s , and their ratio.
We can also estimate m u (or equivalently the ratios m u /m d or m u /m) from our simulation.
Given m s , we find m u by extrapolating in the light valence mass to the point where the K + has the mass (m K + ) QCD , where "QCD" indicates that EM effects have been removed. We with δ E = 1 our central value, corresponding to ∆ E = 1 and vanishing EM correction to the  . IX C and Table IV), we find a negligible error ≈ 0.01%.
Our simulation directly determines decay constants of the charged mesons, π + and K + , in the absence of electromagnetism and with m u = m d . Since mπ is (approximately) the π + mass in this limit, we must simply extrapolate f π + to the point where the mesons have the masses in Eq. (45), i.e., to our physical values ofm and m s . The situation with the kaons is rather different. It is f + K that is measured experimentally, not some isospin averaged decay constant of K + and K 0 . We therefore should extract f + K by extrapolating the light valence quark to the physical value of m u , notm. Despite the large uncertainty in m u /m d from EM effects, the indirect error induced in f + K through m u is tiny, ∼ 0.07%. This is due to the fact that the decay constant changes slowly with mass: It varies only by 20% all the way from the K to the π.
In principle there are also direct isospin-violating errors in the decay constants. For f K + , there is an effect we can estimate from the coefficient β

IX. S χ PT FITS
We fit the partially quenched data for m 2 P + 5 /(m x + m y ) and f P + 5 together in all fits; this helps to constrain the common O(a 2 ) chiral parameters. Similarly, both coarse and fine data are fit together, helping to constrain the overall lattice spacing dependence. Correlations between and among masses and decay constants within each sea-quark set are included, with the covariance matrix computed as described in Sec. III.

A. Data subsets and fit ranges
Our lattice data is very precise (0.1% to 0.7% on m 2 P + 5 /(m x + m y ), and 0.1% to 0.4% on f P + 5 ); while the χ PT expansion parameter for the kaon, χ ud,s , is ≈ 0.18. Since χ 2 ud,s ≈ 0.03, we cannot expect NLO χ PT, which is missing corrections to m 2 P + 5 /(m x + m y ) or f P + 5 of order χ 2 q , to work well for meson masses that are even an appreciable fraction of m K . NNLO χ PT, however, may allow us to fit up to fairly near m K , because the missing corrections at the kaon, χ 3 ud,s ≈ 0.006), are comparable to (but somewhat larger than) the statistical accuracy of our data. Of course, this is only a rough guide. There are at least two sources of complications: (1) It is an idealization to imagine that the chiral expansion is governed by a single mass parameter. Many different mesons contribute to chiral loops. Although we can restrict the valence masses in the fit, the s sea quark mass in the simulations is fixed at m ′ s . Thus there will always be some contributions from fairly heavy mesons. (2) Taste violations produce additional contributions to meson masses, or, effectively, add another expansion parameter, χ a 2 , Eq. (12).
In practice we consider three different subsets of our complete (coarse and fine) partially quenched data set. Compared to the strange sea quark mass in the simulations, m ′ s , we can tolerate somewhat heavier valence masses on the fine lattices, since on those lattices m ′ s exceeds m s by a smaller amount and contributions to meson masses from taste splittings are smaller. The sets are: • Subset I . m x + m y ≤ 0.40m ′ s (coarse), and m x + m y ≤ 0.54m ′ s (fine). 47 valence mass combinations; 94 data points.
There are always twice as many data points as mass combinations because we are fitting m 2 P + 5 /(m x + m y ) and f P + 5 simultaneously.
On the fine lattices, m ′ s is about 10% larger than m s , so we expect errors of order (0.54 × 1.1×0.18) 2 ≈ 1.1% at NLO in subset I , without even considering the effects of mixed valencesea mesons or taste violations. Thus we do not anticipate that the NLO form, Eqs. (17) and (18), can fit the data, even on subset I . Good fits should require at least the NNLO forms, Eqs. (32) and (33), on all sets. Indeed, fitting with Eqs. (17) and (18)  We note that it is not practical to use valence masses below those in subset I because we rapidly run out of data, and in any case we cease to reduce significantly the masses of mesons made of a valence quark and a strange sea quark.

B. Inventory of parameters and alternative fits
Since there are a large number of fit parameters, we provide an inventory before discussing the final fits in more detail. Our standard NNLO fits on subsets I and II have the following number and types of parameters:  We remind the reader that the NLO taste-violating parameters [(c) above] are forced to change by a fixed ratio (in a given fit) in going from the coarse to fine lattices. The point is that these parameters start at O(α 2 S a 2 ), so we know how they change with a, up to corrections that are higher order in α S and/or a. A range for the ratio is considered in assessing the systematic error (see discussion following Eq. (10)).
The priors restricting the parameters governing lattice spacing dependence [(f) above] require further explanation. We note first that it is not possible, at least with the current data set, to remove these restrictions on all the physical parameters, allowing them to be arbitrarily different on coarse and fine lattices. If we do that, the fit becomes unstable because there are directions in parameter space in which the fit function is almost flat.
Some of these directions can easily be seen in Eqs. (17) and (18). For example, a 2 L ′ + a 2 L ′′ can grow large, compensated by a decrease in µ and corresponding increases in 2L 8 −L 5 , in Eq. (32). Only the first (continuum) log term in Eq. (17) is uncompensated, but we already know that the good fits allow a fairly large range in its coefficient, the chiral coupling. (See discussion of the parameter ω in Sec. VI B.) A similar mode involves a 2 L ′′ − a 2 L ′ growing and f decreasing in Eq. (18). Even when such modes eventually converge, the resulting fit is completely unphysical, with 20%-100% variation of results with lattice spacing and enormous NLO chiral corrections. Once the physical parameters are required to change by only a small amount with a, however, they cannot compensate for changes in taste-violating parameters like a 2 L ′ and a 2 L ′′ , and the runaway modes are damped.
Since the taste-violating NLO parameters (a 2 L ′ and a 2 L ′′ ) are included explicitly, the LO parameters f and µ should have only generic (∼ 2%) changes with a, which is our preferred choice for priors for their variations. Of course, as mentioned following Eq. (42), differences in µ between coarse and fine lattices can also be due to perturbative errors in the ratio R m . When we use the one-loop value for R m , 0.958, the fits prefer an ≈ 8% difference in the renormalized µ value between coarse and fine, which is one indication of the size of the perturbative error.
For most of the NLO physical parameters [(b)], the preferred 2.5% prior for latticespacing dependence is not restrictive, with the fits finding a change that is significantly smaller ( < ∼ 0.5%). The exception is 2L 8 −L 5 , for which the 2.5% prior results in almost a 6% difference (a 2.3σ effect), suggesting that the corresponding NNLO taste-violating term has a sizable coefficient. If we instead remove any restriction on the a-dependence of 2L 8 −L 5 (while keeping the constraints on the other physical parameters), 2L 8 −L 5 varies by ≈ 20% from coarse to fine, which is sizable, ∼ 3(χ coarse a 2 − χ fine a 2 ). Perhaps generic and taste-violating effects are both contributing significantly in the same direction. Fortunately, the continuum extrapolated value of 2L 8 −L 5 changes by only 5% when we remove the restriction on its a-dependence. In any case, the systematic errors on 2L 8 − L 5 (as well as on the other L i ) are dominated by the larger changes caused by varying the mass range and/or the details of the chiral fits. Thus, even if we were to take the full 20% variation from coarse to fine as the "discretization error" on this parameter the final errors quoted in Sec. X would hardly change at all.
Our preferred 2.5% prior for the NNLO physical parameters [(d)] is again generally not restrictive, with β (m) 2 the only exception. We note that allowing the NNLO terms to vary with a is not systematic, 8 because it effectively introduces some, but not all, NNNLO terms.
Therefore we consider an alternative in which NNLO physical parameters do not vary with a ("0% priors"), but all other features of the preferred fit are unchanged. This fit has lower confidence level (CL = 0.36) than the preferred fit (CL = 0.65), but gives physical results that are very similar.
The alternative choices of priors for parameters (f) also give good fits, except in the cases of 1% (CL ≈ 5 × 10 −3 ) and "0%" priors 9 (CL ≈ 10 −4 ). However we keep all the choices in the systematic error analysis. The 0% priors case, where all physical parameters are fixed as a function of a, gives results that are in no way extreme among all alternative fits considered in this work.
The preferred NNLO fit with parameters (a)-(f) above results in good confidence levels on data subsets I and II . This is all we need to extract the low energy constants L i , and it is acceptable for determining f π andm. But to determine f K and m s without an extrapolation to heavier valence masses, we need to fit to the data in subset III . We can then interpolate to the physical strange mass. However the preferred NNLO fit, or variants thereof, does not give good confidence level when applied to subset III . Based on the discussion in Sec. IX A, this is not surprising, especially considering that the simulation mass m ′ s is larger than the physical strange quark mass. The problem, though, is that it is not possible to move beyond NNLO in anything approaching a systematic way without introducing an unwieldy number of new parameters. If instead we fit subset III to our rather ad hoc NNNLO forms in Sec. VI C, we can obtain acceptable fits. But the high quark masses involved, as well as the nonsystematic treatment of the higher order terms, may introduce significant systematic errors in the low energy constants and therefore in the extrapolation ofm to the physical point, which is, of course, still needed for finding m s and f K . Thus such fits are not acceptable for finding decay constants and quark masses.
Our solution to the above dilemma is to use the results for all LO and NLO parameters from the fits to subset II as inputs to the fits on subset III . We then use the form Eqs. (36) and (37) on subset III , but with the LO and NLO parameters, and their change with lattice spacing, constrained to vary by at most their statistical errors around their previously determined values. The NNLO and NNNLO parameters are left unconstrained. Their variation with lattice spacing, which implicitly introduces higher order mixed terms in χ q and χ a 2 , either is also left constrained or is constrained mildly, with 1-σ constraints of 10% (preferred fit), 15%, or 20%. The total number of parameters in this fit is 48, of which 20 (16 LO and NLO parameters and their a-dependence, plus the usual 4 scale parameters) are 9 In this case there are 24 free parameters in the fit. tightly constrained. We get reasonable fits with all these versions. The advantage of this approach is that we can interpolate to the physical strange quark mass and extrapolate to the physical light quark mass in the same fit. We therefore use the preferred fit to subset III for our central values for quark masses and decay constants. The ad hoc nature of the higher order terms in this fit is not a problem because we have already satisfied the chiral constraints at low quark mass. All that is required in the region of m s is a fit that interpolates from our partially quenched valence masses and nominal sea quark mass to the physical m s value. In other words, we do not need to rely in detail on chiral perturbation theory for the m s dependence, since we can reach the physical value of m s in the simulation.
Effectively, we are depending only on two-flavor χ PT. Note however that we include the more conventional chiral fits on subsets I and II among the alternatives when estimating the systematic errors. Using the volume-dependence from NLO S χ PT, Eqs. (13) and (14), the leading finite volume effects can be removed from our data. Such effects are rather small to begin with (< 0.9% on M 2 P + 5 and < 1.4% on f P + 5 , based on fit B), and this calculated volume-dependence is consistent with simulation results in the one case where two different volumes are available [19]. One-loop finite volume effects have been removed from the points and lines shown in To extract continuum results for masses or decay constants from Fits A, B, or C, we first set the taste splitting and the taste violating parameters to zero. We then extrapolate the remaining, physical parameters to the continuum linearly in α S a 2 . For central values, we assume that the ratio of this quantity between fine and coarse lattices is 0.427 (see Eq. (26)).

C. Fit results
For the LO parameters we take the range of the ratio to be 0.398 to 0.441 in estimating the systematic error, as in the discussion of Eq. (26). But for all other parameters, we expand the range to 0.30 to 0.441 in recognition of the fact that the fits do not distinguish generic discretization errors, O(α S a 2 ), from taste-violating errors, O(α 2 S a 2 ). We thus must include the range discussed following Eq. (10). Table IV shows the central values of the continuum extrapolated parameters for Fits A, B, and C. Note that the statistical errors on most of the parameters are quite large. This seems to be a consequence of the "flat directions" in the fitting function, as described in Sec. IX B: small fluctuations in the data can produce large variations in the parameters. However, because of the correlations among the parameters, the statistical errors of interpolated or extrapolated decay constants and masses are small, comparable to those of the raw data.
Once the continuum chiral parameters are obtained, we set valence and sea quark masses equal and obtain "full QCD" formulas for m 2 π /(2m ′ ), m 2 K /(m ′ + m ′ s ), f π , and f K in terms of arbitrary quark massesm ′ and m ′ s . The cyan lines in Figs. 11 -15 show these as a function ofm ′ , with m ′ s held fixed at the value of the simulation sea quark mass on the fine lattices. As a consistency check, we also show in each case the result from extrapolation of full QCD points to the continuum at fixed quark mass (cyan fancy squares). To generate these points, we use the chiral fits only to interpolate the coarse data so that it corresponds  to the same physical quark masses as the fine data. There are just two such points in each plot because we have just two runs with different sea quark masses on the fine lattice.
Since discretization errors come both from taste violations and generic errors, there is an ambiguity in the extrapolation used to find these points. We have assumed that taste violations dominate and have extrapolated linearly in α 2 a 2 , i.e., with a ratio of 0.35 in α 2 a 2 between coarse and fine (see discussion following Eq. (10)). This agrees both with our order of magnitude estimates (taste violations ≈ 6%; generic errors ≈ 2%) and with a detailed analysis below (Sec. IX D 1).
To proceed further we need to know the physical values of the quark masses. These can be obtained from Fig. 12 or Fig. 15 by finding those values ofm and m s that give the π and K their physical QCD masses in the isospin limit, mπ and mK (defined in Eq. (45)). An iterative procedure is required because both meson masses depend on both quark masses, although the dependence of mπ on m s is mild, since s only appears as a sea quark. The nature of this extrapolation/interpolation is most clearly seen in Fig. 16, 10 which is again Fit We then fix the sea quark mass atm and continue the extrapolation in valence mass m x to m u (short green dotted line). It is clear from the size of the systematic errors that this final short extrapolation does not at present produce a significant change in f K . However, this distinction between f K + and fK will become more important as lattice computations improve.

D. Discussion
Good fits are not possible without the taste-violating terms in S χ PT. Figure 10 corresponds exactly to Fig. 11 except that in Fig. 10 the taste splittings in meson masses have been set to zero and the taste-violating chiral parameters (δ ′ V , δ ′ A , L ′ , and L ′′ ) have been eliminated. Thus the fit in Fig. 10 is to the "continuum" NNLO χ PT form, 12 which has only four fewer parameters than the S χ PT fit in Fig. 11 (i.e., a total of 36). We put "continuum" in quotation marks here because generic variations in physical parameters between coarse and fine lattices are still allowed. Further, even if we allow these generic variations to be 11 See Sec. X for additional discussion about the experimental value of f K . 12 As in the S χ PT case, NNLO chiral logarithms are not included.
arbitrarily large, instead of the ≈ 2% variation permitted in the standard fits, we cannot obtain good fits without including the taste violations. The fact that continuum χ PT fits are so poor reassures us that the good fit in Fig. 11 is not simply a trivial consequence of having a lot of fit parameters -one has to get the physics right. Other test fits described below, such as fitting without the nonanalytic terms in the fit function (Sec. IX D 2), give additional reassurance, since they have equal or comparable numbers of parameters to Fit B but cannot describe the lattice data.

Continuum extrapolation
For the decay constants, our preferred method of continuum extrapolation is to extrapolate the chiral fit parameters, as described in Sec. IX C. An alternative method is to determine decay constants in physical units at fixed lattice spacing, and then attempt to extrapolate these quantities. There are two ways to find fixed lattice-spacing values: (1) simply use the complete chiral fits to extract the decay constants on the coarse or fine sets, or (2) first set the taste-violating parameters (δ ′ V , δ ′ A , L ′ , L ′′ ) and splittings in the fit to zero, and then extract the decay constants for each set. The advantage of method (2) is that, once taste-violations have been set to zero, remaining discretization errors should be dominantly of the generic type, so we may extrapolate to the continuum linearly in α S a 2 . In method (1), the decay constants at fixed a have both generic and taste-violating discretization errors, so there is an ambiguity in extrapolating them to the continuum. Figure 18 shows the decay constants at fixed a using both methods, and the various extrapolations to the continuum. Once taste violations are removed (method (2)), the remaining discretization errors are quite small, giving a 1%-2% change between coarse and fine, as expected. Further, because this change is small, it would not make a significant difference in the extrapolated answer if we were to replace the α S a 2 extrapolation with an α 2 S a 2 one. With method (1), the change from coarse to fine is 5%-6%, roughly the same size as the difference in the raw data between these lattices: See Figs. 11 -15. 13 There is also a noticeable difference (≈ 1.2%) in the extrapolated results with method (1), depending on whether O(α S a 2 ) or O(α 2 S a 2 ) errors are assumed. This ambiguity would grow to ≈ 2.2% if we used the full allowed range of values for the O(α S a 2 ) or O(α 2 S a 2 ). For f π and f K , we therefore do not consider method (1) results among the possible alternatives in assessing the systematic error. We use the parameter-extrapolated version (red diamond and fancy plus in Fig. 18) for central values, and include the results of method (2) when estimating systematic errors. Note that we are rejecting method (1) because of its inherent ambiguity, not because it disagrees with the other methods of extrapolation.
On the other hand, for the ratio f K /f π and for quark masses (m u /m d here, and m s ,m and m s /m in Ref. [22]), the results change little with lattice spacing, so the ambiguity in method (1) is tiny (much less than other systematic errors). Therefore, in those cases we include all three methods of extrapolation in our systematic error estimates.
The nice consistency of Fig. 18 with our understanding of the sources of discretization errors is comforting. However, we caution the reader that some aspects of this picture are dependent on the assumptions that went into our fits. In particular, we have inserted Bayesian priors to enforce a (1 standard deviation) constraint that LO chiral parameters change by at most 2% from coarse to fine lattices. When we relax this constraint to 4%, the difference of method (2) results from coarse to fine increases to 3%-4%; while method (1) differences remain at about 6%. The relaxed version of the chiral fit is included in the systematic error analysis. With just two lattice spacings, however, we cannot remove this constraint entirely without the fit becoming unstable, as has been emphasized in Sec. IX B.
Instead, the key point here is that we can obtain good fits of our entire data set based on the theory of taste violations, plus some smaller generic errors. This is to be contrasted with our failed attempts to fit the data without including taste-violations, even when generic errors are allowed to be arbitrarily large.

Evidence for chiral logarithms
From Sec. VI B, we know that the coefficient of the chiral logarithm terms, the chiral coupling 1/(16π 2 f 2 ), is not tightly constrained by the fits. If the chiral coupling is allowed to be a free parameter, fits to higher valence masses prefer values of f near f K in the coupling; while fits to lower masses prefer f closer to f π . On the other hand, acceptable fits can be obtained for all our mass ranges with the chiral coupling fixed anywhere between its value for f π and that for f K . Given this freedom, can we claim that chiral logarithms are observed at all? To answer this question, we consider a variety of alternative fits without chiral logarithms. Bećirević and Villadoro [46], have pointed out that, for some current simulations on small volumes, the finite size effects are much more important than the actual chiral logarithms.
We would not expect that to be the case here since the lattice volumes are relatively large (L ≥ 2.5 fm) and finite volume effects here are small (at most 1.4% -see Sec. IX C). To check this expectation, we first removed the finite volume effects from the data using Fit B, and then fit again to the 38 parameter form with ℓ =l = 0. This fit is improved over the previous one, but still quite bad: χ 2 /d.o.f. = 3.08 for 202 degrees of freedom; CL < 10 −43 .
Our conclusion is that the effect of the staggered chiral logarithms is in fact observed in our data.
One can ask whether the finite volume effects are also directly observed. The answer seems to be yes: A 40-parameter fit leaving out these effects (setting δ 1 and δ 3 to zero in Eqs. (13) and (14)) but otherwise identical to Fit B has χ 2 /d.o.f. = 1.95 for 200 degrees of freedom; CL ≈ 2 × 10 −14 .

Are the lattice masses light enough for S χ PT to be applicable?
To discuss this question, we first have to say what we mean by the physical quark masses at fixed a. For current purposes we define the physical values of the lattice masses, am s and am, by method (2), i.e., as the quark masses that give the pion and kaon their physical masses when all taste splittings and taste-violating chiral parameters are set to zero. This gives am s ≈ 0.0457, am ≈ 0.00166 on the coarse lattices and am s ≈ 0.0289, am ≈ 0.00105 on the fine. We have not made a detailed study of the errors in these numbers, but systematic errors are ∼ 6%, and statistical errors are 1% or less. We could alternatively define the physical m s at a given lattice spacing by method (1), i.e., as the masses that give the physical mesons masses directly, including all effects of taste violations in the chiral loops. The latter values of quark masses were quoted in [19,22] and are about ∼ 15% smaller on the coarse lattices and ∼ 6% smaller on the fine lattices than the ones quoted above. We choose the method (2) definition here because we are going to be adding on the taste-splittings to meson masses explicitly. We note that the difference between the methods becomes becomes small, 1%-2%, once we extrapolate to the continuum, and is included in the systematic errors quoted in Ref. [22].
We now consider meson masses on the coarse lattice for subset II . This is the worst case because subset II contains the highest valence quark masses for which we have applied a fully chiral description, and the coarse lattices have the largest additional contributions to meson masses from taste violations. The valence masses here obey m x + m y ≤ 0.7m ′ s ≈ 0.77 m s . Since the smallest valence mass is 0.1 m ′ s , the largest is 0.6 m ′ s ≈ 0.66 m s ; while the largest sea quark mass in the simulation is m ′ s ≈ 1.09 m s . For the following estimates, we assume linear dependence of squared meson masses on quark masses, and take 486 MeV as the mass a "kaon" would have in the absence of electromagnetism and if the light quark were massless. In other words, we use µm s = (486 MeV) 2 . This comes from Eq. (45) and the ratio m s /m = 27.4 [22]. Then the largest valence-valence Goldstone meson mass in coarse subset II is 425 MeV. Adding on the largest splitting (the taste-singlet case) gives 623 MeV; while an "average" taste splitting (see discussion after Eq. (12)) gives 551 MeV. We do not think it unreasonable to expect S χ PT to work in this mass range, although it is not surprisingconsidering the small statistical errors of the data -that NNLO terms are needed. Further, the comparison with subset I results is a good check, because the corresponding masses there are significantly lower: 321 MeV, 557 MeV, and 475 MeV, respectively.
Mesons with one or two sea quarks also appear in chiral loops. These are generally comparable to the masses just discussed. But they are significantly larger when the sea quark is an s, which is exacerbated by the fact that the simulation value, m ′ s , is larger than the physical mass. On subset II , this largest valence-sea (Goldstone) mass is 642 MeV.
Adding on the biggest splitting gives 787 MeV. This taste-singlet meson enters with a factor 1/16 in the sum over tastes; the "average" taste version is 731 MeV. On subset I , the masses become 578 MeV, 736 MeV, and 676 MeV, respectively. Further, there are sea-sea contributions, which are independent of the valence mass subset. The most relevant here is the taste-singlet η. Its mass depends (mildly) on the light quark sea masses, but is ∼ 765 MeV including splitting. In addition there are η ′ -like particles in the taste axial and vector channels whose masses are comparable to the taste-singlet η but have smaller (O(a 2 )) couplings.
The meson masses involving the s sea quark are admittedly quite high to expect that even NNLO χ PT will be accurate. For example, the largest mass mentioned above, 765 MeV, corresponds to a χ q value, Eq. (11), of 0.43. This suggests an error from neglected terms of order (0.43) 3 = 8%. But, just as for the valence subset III, the issue here for decay constants and quark masses is only to get a good interpolation to the physical s quark mass. Indeed, if the s sea quark in the simulation had been chosen at the a posteriori determined physical mass, we would not have needed to use χ PT for the s at all, but could use a SU L (2)×SU R (2) χ PT for the light quark extrapolation. The systematic error on the coarse lattice from adjusting m ′ s = 1.09m s to m s may be crudely estimated as (0.43) 3 × [(1.09) 3 − 1] ≈ 2%. Since m ′ s is closer to m s on the fine lattice, some of this error will be extrapolated away when we go to the continuum limit. On the other hand, chiral coefficients (at NLO and NNLO) that involve the sea quarks are not fit accurately because the "lever arm" is small: the sum of the sea quark masses changes by less than a factor of 2 over our entire range of coarse lattices and only by 30% for the fine lattices. A more reliable estimate of the error in adjusting the s quark mass comes from considering the range in results over the full list of alternative mass subsets, chiral fits, and continuum extrapolations. It can be as large as half the total chiral error in our results for decay constants and quark masses (see line a1 in Tables V and VI); the remaining error comes from extrapolating in the light quark mass.
For the L i , the situation is somewhat different. Missing higher order terms in the SU L (3) × SU R (3) expansion mean that there is spurious analytic dependence on the light quark masses that increases as meson masses get larger. Here we are missing the NNLO chiral logarithms, so those terms determine the size of the errors. Letting M be a generic meson mass, the absent terms are of order M 4 log(M 2 )/(8π 2 f 2 π ) 2 . Putting in the largest meson masses discussed above, results in an estimate of the absolute error in the L i of a few times 10 −4 . This is indeed the size of the errors we observe when we consider all the alternative chiral fits discussed above and/or restrict the valence masses to subset I instead of subset II (Sec. X). The low energy constant 2L 6 − L 4 may be an exception: Since the errors in it are large on this scale, ∼ 4 × 10 −4 , they may also not be very reliable. Difficulty in extracting 2L 6 − L 4 is again related to the small lever arm for the sea quark dependence.
This can be seen in Fig. 19, an enlargement of a small region of Fig. 12. As the sea quark mass is changed, the differences are small -comparable to statistical errors -and not monotonic. Contrast this with the monotonic sea quark dependence seen for f π in Fig. 11.
A coarse simulation now in progress, with all three sea quark masses at about 0.66m s , should help to reduce significantly the error in the sea quark mass dependence. Figure 20 shows the convergence of SU(3) L × SU(3) R χ PT for f π and f K . All chiral parameters in this plot have been extrapolated to the continuum. The NLO terms contribute 20%. This is true even for f π , because m s does not vanish in the chiral limit. The convergence of SU(2) L × SU(2) R χ PT for f π can also be extracted from this plot by starting with the "chiral limit" line instead of the "LO" line lowest order contribution to f π . Note that the magenta and cyan fancy squares, which are included as a consistency check, are the only full QCD points that we can extrapolate to the continuum at fixed mass. We have lighter valence quarks on the fine lattices, and lighter valence and sea quarks on the coarse lattices.

Convergence of χ PT
All such partially quenched points are included in the fits that produce the lines in the plot.  Table IV), so we must accept the large systematic effect on 2L 6 − L 4 as inherent in the current data set.
Indeed, the size of the difference in the LO term, µ, between the two fits is reasonable, given that the fits prefer a chiral coupling 1/(16π 2 f 2 ) with f moving from ≈ f π to ≈ f K as the quark masses rise toward m s (see comments about the parameter ω in Sec. VI B). Both fits here fix f = f π , but Fit B is effectively able to reduce the effect of the chiral coupling by reducing the LO parameter µ and compensating by increasing the NLO parameter 2L 6 − L 4 .
Since the difference between f 2 π and f 2 K is more than 40%, the ∼ 20% difference between the µ from the two fits is not unexpected. New simulations with lighter strange sea quark masses will allow us to take all quark masses deeper into the chiral regime, as well as greatly increase the lever arm on the m s dependence, and should help to resolve this issue.

Mass dependence of renormalization scheme
Strictly speaking, it is incorrect to use a scale determined by a quantity like r 1 or the Υ ′ -Υ mass difference in chiral fits to lattice data, since such quantities themselves have some (small) sea quark mass dependence not included in chiral perturbation theory [47]. We have investigated this effect by changing to a mass-independent renormalization scheme: Instead 14 The statistical errors here are slightly larger than those in Table IV because the statistical errors associated with the continuum extrapolation are included. of fixing the (relative) lattice scale for particular sea quark mass values from r 1 /a at those values, we can use, on each lattice, the value of r 1 /a after extrapolation to the physical mass values 15 using the fit Eq. (7). This produces the same value of a on all the coarse (or, separately, fine) lattices, independent of the sea quark masses. We include the difference between the mass independent scheme and our standard approach in our assessment of systematic errors.
For decay constants and quark masses, the change in scheme is a priori unlikely to make much difference because the physical point is unaffected -all that may change is the extrapolation to it. Further, any low-order analytic dependence on sea quark mass introduced though r 1 would automatically be compensated by changes in the analytic chiral fit parameters. So the only problem would be due to nonanalytic quark mass dependence, which is probably quite small because r 1 is a short-distance quantity, at or near the perturbative region. We thus consider the variation in scheme simply as another alternative version of the chiral fits. This means that it would affect the final systematic error only if it produced the largest difference from the central value over all the alternatives. In fact it is fairly small, as expected (see line a3 in Tables V and VI).
The situation is logically quite different for the low energy constants. Here, analytic sea quark mass dependence in r 1 would directly change the output values of L 4 and 2L 6 − L 4 , which multiply sea quark masses. We therefore consider the scheme dependence as a systematic error in its own right, and add any error found in quadrature with other systematic errors. In practice, however, this effect is still smaller than other errors (see Table VII).

Residual finite volume effects
At the precision we are working (especially for f K /f π ), it is important to consider whether finite volume effects coming from terms beyond one-loop in S χ PT could be non-negligible.
Indeed, Colangelo and Haefeli [48] have recently investigated such effects in full continuum QCD. For volumes and masses comparable to those used here, they find large higher order corrections to the volume dependence, roughly 30% to 50% of the one-loop results. 15 Here we consider both the quark mass values determined by method (1) and those determined by method (2) -see Sec. IX D 3.
In asymptotically large volumes, the finite-volume effects in S χ PT are suppressed relative to those in continuum χ PT for the same (Goldstone) masses because most of the pions entering chiral loops have larger masses. This can be easily seen in Eq. (24): the lightest (Goldstone) pion appears with a weight 1/16 as large as in the continuum. However, at our current volumes, masses, and lattice spacings, the relation between S χ PT and χ PT finite volume effects is complicated, with the former just as likely to be larger than the latter as smaller. Equations (24) and (25), together with our numerical results for splittings (Table III) and the taste-violating hairpin parameters (Eq. (47), below), show how the asymptotic rule can be violated. Since ∆ A is the smallest splitting, and δ ′ A is non-negligible and negative, m η ′ A may not be much larger than the Goldstone pion mass. Then, due to the factors of 4 in Eq. (24), finite volume corrections coming from the η ′ A can be as large or larger than the continuum corrections. As a → 0, the term ℓ(m 2 η ′ A ) would be canceled by ℓ(m 2 π A ), but this cancellation may not be effective for finite volume effects at a given value of a, because the volume effects are sensitive to small mass differences. Since one-loop finite volume effects on our lattices are comparable to those in the continuum, we have no a priori reason to expect that our results are protected from the higher order effects [48] seen in the continuum.
Assuming that the higher order finite volume effects in S χ PT are roughly the same size as those in full continuum QCD when the volumes and the masses of the mesons in the loops are the same, we can estimate the resulting systematic errors. For our data, the biggest one-loop finite volume effects appear when both the valence masses andm ′ are small (giving a light η ′ A in the two-flavor case, Eq. (25), or a light η A in the three flavor case, Eq. (22)). The worst case occurs in the coarse run with am ′ = 0.007; the run withm ′ = 0.005 has smaller finite volume effects because L ≈ 3.0 fm there, instead of L ≈ 2.5 fm for other runs.
From the calculations in Ref. [48], 16 we estimate that the residual higher order finite volume effect is at most 0.47% in f π , 0.24% in f K , and 0.23% in f K /f π . More stringent bounds on the errors can be obtained by removing from the data set those points that have the largest finite volume corrections. Eliminating 8 of 240 points from mass subset II (5 from coarse run 0.007/0.05, 2 from coarse run 0.01/0.05, and 1 from fine run 0.0062/0.031), we lower the 16 We are indebted to Gilberto Colangelo for providing us with the results for higher order effects in f π and f K at the values of volume and meson mass relevant to our computations largest one-loop finite volume effect on f P + 5 or M 2 P + 5 from 1.35% to 0.81%. Not surprisingly, since the reduced data set retains most of the lowest valence mass points and all of the lowest sea mass points, it produces nearly identical results (within 0.05%) as the original set. But a repeat of the analysis using Ref. [48] now bounds the residual finite volume error by 0.29% in f π , 0.15% in f K , and 0.14% in f K /f π . These are negligible compared to our other systematic errors. In the future, however, as quark masses in staggered simulations decrease further, it will be necessary either to have a better handle on these higher order effects in S χ PT or to go to significantly larger volumes.

Fourth root of the determinant
In order to eliminate the quark doubling that is still present in the staggered action, the simulations here take the fourth root of the quark determinant for each flavor in order to reduce the quark tastes from 4 to 1 per flavor. There is apparently no ultra-local lattice action that would correspond to the effective action that results from taking this fourth root.
The possibility thus exists that physical non-localities will remain in the continuum limit, potentially spoiling the description of QCD by the staggered action. The good agreement of the staggered results with experiment and with continuum chiral behavior plus understood discretization effects (both in current and previous work [17,19,20]) lead us to believe that this is not a problem, but the question is not settled.
The comparison of simulation data with S χ PT forms allows us to make a crude but somewhat more direct test of the fourth root trick. Equations (17) and (18), as written, take into account the fourth root by dividing each sea quark loop contribution by 4, to leave 1 taste per flavor. It is a simple exercise to generalize Eqs. (17) and (18)  on subsets II and I , respectively. In the latter case, the coefficient ω corresponds to an f in the chiral coupling that is about halfway between f π and f K (ω ≈ 0.82(11)), which is reasonable. Given that the fits in any case do not tightly constrain the chiral logarithm terms (see Sec. VI B), we consider these results satisfactory.
We note that there has been recent numerical [49] and analytic [50] work indicating even more directly that the fourth root trick should work as expected. On the other hand, there have been two other recent papers that purport to show problems with locality [51].
We do not believe the latter work is worrisome because it does not take into account the taste structure of staggered quarks. Instead of trying to project onto a single taste to find the fourth root of the determinant, those papers look only at the fourth root of the Dirac operator itself. That procedure, in our opinion, is almost guaranteed to find a nonlocal result, just as it would in trying to reduce eight Wilson fermions to two, which certainly has an alternative, local solution.

Taste violating hairpins
Before turning to our physical results, we quote the values of the two taste-violating hairpin parameters coming from the fits. Together with the splittings, Table III, these parameters appear in S χ PT calculations for other physical quantities, such as heavy-light decay constants [52]. Averaging values from Fits A and B, we find, on the coarse lattices: where the errors are statistical and systematic, respectively. The latter error comes from the variation over all acceptable chiral fits on mass subsets I and II . The parameter a 2 δ ′ A is comparable in size to the taste-violating splittings (Table III); while a 2 δ ′ V is consistent with zero but poorly determined. The values of a 2 δ ′ A and a 2 δ ′ V on the fine lattices are not fit separately but are constrained to be 0.35 times as large for central-value fits. (See discussion following Eq. (10).)

X. FINAL RESULTS AND CONCLUSIONS
The central values and error estimates for f π , f K , and f K /f π are collected in Table V.
Central values come from Fit C (mass subset III ). The scale errors are found by repeating the analysis after moving our value r 1 = 0.317 fm by plus or minus one standard deviation (±0.007 fm); see Sec. III. The change in decay constants under this variation in the scale is slightly less than the nominal 2.2% = 0.007/0.317. This is due to a cancellation coming from the corresponding readjustment of the quark masses needed to give the mesons their physical masses.
The (indirect) EM errors just come from changes in our results for quark masses due to the assumed range of ∆ E (see Sec. VIII); clearly this effect is very small. Direct EM effects that pertain to the comparison of decay constants with experiment are much larger, of order several percent -see Ref. [37]. However the direct effects are not relevant here because f π and f K are defined in the absence of electromagnetism.
The chiral/continuum errors are found by taking the maximum deviation from the central value over all versions of the chiral fits described in Sec. IX B and Sec. IX D 5, and all versions of the continuum extrapolations described in Sec. IX C and Sec. IX D 1 (including ranges in assumptions about how α S a 2 and α 2 s a 2 change from coarse to fine lattices -see discussions Eqs. (26) and (10)), as well as variation in the perturbative parameter R m described in Sec. VII. Because the continuum and chiral extrapolations are connected within S χ PT, it is not meaningful to quote separate errors for each. However, since a large number of alternatives are considered here, we believe it will be helpful to the reader to report the variations in physical results as one moves along various "slices" through the alternatives.
The slices shown in Table V are defined as follows: a. All alternative chiral fits on all mass subsets, but only with the preferred method of continuum extrapolation (extrapolation of chiral parameters), and only with preferred values of the ratios of α S a 2 and of α 2 S a 2 (Eqs. (26) and (10)).
a1. Same as a, but restricted to mass subset III. This is mainly an estimate of the errors involved in interpolating around m s .
a3. Same as a, but restricted to fits where the scale is chosen in a mass-independent manner (see Sec. IX D 5).
b. Alternative values of the ratios of α S a 2 and of α 2 S a 2 used in continuum extrapolation and/or alternative method of extrapolation (method (2) -Sec. IX D 1) and/or  alternative value of R m (Sec. VII). The preferred chiral fit is kept (Fit C).
b1. Same as b, but restricted to the preferred value (0.35) of ratio of α 2 S a 2 for all taste-violating quantities.
b2. Same as b, but only the ratio of α 2 S a 2 is varied (in the range 0.3-0.4) and only for taste-violating quantities that are not directly measured (δ ′ A , δ ′ V , L ′ , and L ′′see Sec. VI A). The preferred continuum extrapolation (extrapolation of chiral fit parameters) is used.
b3. Same as b, but only R m is varied, and the preferred continuum extrapolation is used.
c. Alternative method of extrapolation (1) is used and ratio of a 2 varies over union of ranges of α S a 2 and of α 2 S a 2 -Sec. IX D 1.
As discussed in Sec. IX D 1, method (1) continuum extrapolation (slice c) is not included among our systematic alternatives because of the large ambiguity in how to perform the extrapolation. Table V shows, however, that it produces deviations comparable to the full chiral/continuum extrapolation error.
We add in quadrature the signed errors from the chiral/continuum extrapolation, the scale determination, and from direct EM effects, giving a total positive and a total negative systematic error. We then take the larger of the two as a final symmetric error. Note that chiral extrapolation errors and scale errors contribute almost equally to the systematic error on f π and f K ; while scale errors are unimportant for the ratio. The final results for decay constants are: f π = 129.5 ± 0.9 ± 3.5 MeV where the first error is statistical and the second is systematic.
In Sec. VI B we argued that fits that allowed the chiral coupling to vary by more than 10% ("choice (3)" fits with arbitrary ω) should be excluded from the analysis. If we were to include all choice (3) fits in the systematic error analysis, the error on f π would increase from 3.5 MeV to 4.3 MeV; that on f K /f π would increase from 0.013 to 0.022; while that on f K would be unchanged.
Our results are in good agreement with the experimental numbers [37]: f π = 130.7 ± 0.4 MeV, f K = 159.8 ± 1.5 MeV, f K /f π = 1.223 (12). Note that the experimental determination of f K has a rather large error. That is because it depends not only on the precisely measured leptonic decay width of the kaon, but also on V us , which has a significant uncertainty. The errors on our result for f K /f π are small enough that one may turn the comparison around, and use our answer together with the measured leptonic decay widths to constrain V us [2]. With Eq. (16) in Ref. [2], |V ud | = 0.9740 (5), and the current result for f K /f π , we obtain |V us | = 0.2219 (26) .
The error is completely dominated by current lattice errors, which we have added in quadrature. Neglecting |V ub | 2 , the unitarity relation is then The 2σ violation that comes from using the PDG value |V us | = 0.2196 (26) [37] becomes a 1.4σ effect here. We note also that our result is compatible with the very recent KTeV determination [53]: |V us | = 0.2252(8) (21).
The values for f π and f K in Eq. (48) should be considered as updates of those presented in Ref. [17]. The current results are based on an expanded data set. In addition, the analysis in Ref. [17] was performed differently: The data was first extrapolated to the continuum at fixed quark mass and then fit to continuum χ PT forms. S χ PT was used only in estimating the systematic error of the extrapolation procedure. A correction for finite volume effects could not be made with the older approach; instead a finite volume error had to be included.
The present results and those in Ref. [17] agree within their respective systematic errors.
Errors for our direct determination of m u /m d are shown in Table VI. Adding the scale and chiral/continuum extrapolation errors in quadrature, and symmetrizing as for the decay constants, we get the total simulation error. Our final result is where the errors are from statistics, simulation systematics, and direct EM effects, respectively. We have allowed for EM effects in a wide range 0 ≤ ∆ E = δ E ≤ 2 (see Eqs. (45)   Values for quark masses at scale 2 GeV, as well as the ratio m s /m, were reported in Ref. [22]. Since that work used the same lattice data, chiral fits, and error analysis as that described above, we repeat the results here for completeness: where the errors are from statistics, simulation, perturbation theory, and electromagnetic effects, respectively.
Combining the current result for m u /m d with the perturbative mass renormalization calculated in [22] (or, equivalently, withm MS in Eq. (51)), we obtain: where the errors have the same meaning as in Eq. (51), and the scale is again 2 GeV. The separate EM errors in m u and m d are highly, and negatively, correlated, and therefore consistent with the large EM error in m u /m d .
The results for m u /m d and m s /m in Eqs. (50) and (51) appear inconsistent with the relation between m s /m d and m u /m d shown in Fig. 1 of Ref. [4]. However, that appears to be due to NNLO effects not included in [4]. Indeed, Amoros et al. [45] obtain m u /m d = 0.46(9) with a NNLO phenomenological analysis. Further, our results for the two ratios are consistent with the NNLO relation shown in Fig. 3 of Ref. [45].
Since m u is bounded well away from 0, the issue of the physicality of m u = 0 [6] does not arise directly here. However, should the existence of non-perturbative, additive shifts in masses proposed in Ref. [6] be confirmed, there could be some lattice scheme dependence in the quark masses and ratios in Eqs. (51) and (52). We would expect that such nonperturbative effects at the scale of the cutoff would be small at the mass values found here, but there is no proof of this. Comparison with three-flavor results with other lattice regularizations will be important in resolving this question. and IX D 5, we include two additional systematic errors here, to be added in quadrature with the scale and chiral/continuum extrapolation errors: the NNLO error caused by taking µ → µ tree in the NLO terms, and the effect of using a slightly mass-dependent renormalization scheme.
The chiral/continuum "error slices" in Table VII have the same meaning as for the decay constants, except that a1 and c no longer apply. (Slice a1 shows differences with mass set III , which is not included in this part of the analysis, and slice c is not relevant since these quantities are themselves fit parameters.) Further, slice a3, the effects of the mass-dependent scheme, has now been promoted to a separate error.
After adding the systematic errors in quadrature and symmetrizing as before, we obtain: Systematic errors here are dominated by differences over acceptable fits. The chiral scale is taken as Λ χ = m η throughout. Including all choice (3) fits in the systematic error analysis would not change the errors.
The result for 2L 8 − L 5 is well outside the range that would allow for m u = 0 [4,5,7] in the context of χ PT: We note, however, that the constraint on m u coming from 2L 8 − L 5 is not independent give us some control over NNLO effects, the direct determination seems preferable, and can become quite precise if one uses more information on EM effects. This information may come from phenomenology, e.g., Ref. [42], or from lattice simulations, perhaps along the lines of Refs. [55] or [56].
Our approach to computing low energy constants has much in common with earlier work by Nelson, Fleming, and Kilcup [7], who also performed a partially quenched analysis using 3 flavors of dynamical staggered quarks. The main advances in the current analysis are: (1) use of the improved dynamical staggered action and finer lattice spacings, putting us closer to continuum physics, and (2) use of S χ PT to control lattice artifacts, which are still quite large, despite (1). Our result for 2L 8 − L 5 is marginally consistent with that by Nelson et al. [7]; converting their result to chiral scale m η , we get 2L 8 − L 5 = −0.57(1)(14) × 10 −3 .
The current work will be improved by additional simulations now in progress, including coarse lattices at lower strange quark mass (am ′ s = 0.03) and fine lattices at lower light quark mass (am ′ = 0.1am ′ s = .0031). These simulations should enhance our control of the chiral extrapolation, the interpolation around the s quark mass, and the extraction of low energy constants. In addition, we are beginning a parallel analysis on a large quenched data set. If the corresponding S χ PT forms can describe that data well, it will increase our confidence that the interaction of discretization and chiral effects is understood. Beyond that, planned simulations at still finer lattice spacings will provide a better handle on both generic and taste-violating discretization errors, thereby significantly reducing the final systematic errors.  Table I). The red crosses are random-wall source and Coulomb-wall sink, and the green octagons are Coulomb-wall source and point sink (summed over spatial sites to project out the zero momentum states). The blue bursts are from a random-wall source and point sink, and the magenta squares have a Coulomb-wall source and sink. The lower set of "WW" points include an excited state in the fit. The symbol size is proportional to the confidence level of the fit, with the symbol size in the labels corresponding to 50%.
FIG. 2: Same as Fig. 1 but for pion propagator amplitudes. The lower set of "WW" points again include an excited state in the fit. The "PW" symbols have been displaced slightly to the right to separate them from the "WP" points.  Table I).
The amplitudes have been multiplied by 175. Points and fit lines have been corrected for finite volume effects.
Because taste splittings are smaller for the fine lattices, the average meson mass changes more rapidly with quark mass, and there is greater curvature at small quark mass. The magenta squares and pluses have the taste-violating effects at fixed lattice spacing removed with S χ PT ("method (2)"); the points are then extrapolated to the continuum linearly in αa 2 .
The red diamond and fancy plus are the results of extrapolating the chiral fit parameters to the continuum. Extrapolated points at a = 0 have been moved slightly horizontally for clarity. Fit C is used everywhere.