Light Meson Physics from Maximally Twisted Mass Lattice QCD

We present a comprehensive investigation of light meson physics using maximally twisted mass fermions for two mass-degenerate quark flavours. By employing four values of the lattice spacing, spatial lattice extents ranging from 2.0 fm to 2.5 fm and pseudo scalar masses in the range 280 MeV to 650 MeV we control the major systematic effects of our calculation. This enables us to confront our data with chiral perturbation theory and extract low energy constants of the effective chiral Lagrangian and derived quantities, such as the light quark mass, with high precision.


Introduction
In recent years, the non-perturbative description of QCD on the lattice has made significant progress in tackling systematic effects present in the determination of several important physical quantities (see e.g. [1,2] for recent reviews). Simulations containing the dynamics of the light-quark flavours in the sea, as well as those due to the strange quark and recently also to the charm, using pseudo scalar masses below 300 MeV, spatial lattice extents L ≥ 2 fm and lattice spacings smaller than 0.1 fm are presently being performed by several lattice groups. Such simulations will eventually allow for an extrapolation of the lattice data to the continuum limit and to the physical point while keeping also the finite volume effects under control.
Adding a twisted mass term [3] to the Wilson-Dirac operator and tuning the theory to maximal twist [4,5] has by now proved to be a practical and successful tool for performing O(a)-improved lattice QCD simulations, see e.g. refs. [6][7][8][9][10][11][12][13][14][15][16][17]. A detailed study of the continuum-limit scaling in the quenched approximation [18][19][20][21] revealed not only that with an appropriate tuning procedure to maximal twist lattice artefacts follow the expected [4] O(a 2 ) scaling behaviour, but also that the remaining O(a 2 ) effects are small, in agreement with the conclusions drawn in ref. [22]. The only exception seen so far is the neutral pseudo scalar mass [23] which shows significant O(a 2 ) effects. We will discuss this issue and its interpretation [24,25] further in sec. 3.4.1 below.
In this paper, we shall present a study of the continuum limit scaling of lattice QCD with maximally twisted mass fermions for the case of dynamical fermions with N f = 2 mass degenerate quarks demonstrating that also in this setup the O(a 2 ) effects are small as has already been discussed in refs. [26][27][28]. We use data for the charged pseudo scalar meson mass m PS and the decay constant f PS computed in a range of 280 MeV m PS 650 MeV. Our analysis is concentrated on two values of the inverse bare coupling β, β = 3.9 and β = 4.05 corresponding to values of the lattice spacing of a = 0.079 fm and a = 0.063 fm, respectively. We complement our data set by adding an additional coarser lattice spacing of a = 0.100 fm corresponding to β = 3.8 and a finer lattice spacing of a = 0.051 fm, corresponding to β = 4.2. Our new results at β = 4.2, computed at two values of the pseudo scalar mass, confirm the smallness of lattice artefacts in the quantities considered here.
The smallness of lattice artefacts, the range of masses covered, the usage of several lattice volumes and the precision with which m PS and f PS can be obtained using maximally twisted mass fermions allow us to confront the numerical data with chiral perturbation theory and to eventually extract accurate values for a number of low energy constants, in particularl 3 = 3.50 (31) and ℓ 4 = 4.66 (33), as will be discussed below. The main physical results we obtain from this analysis are the light quark mass m MS u,d (µ = 2 GeV) = 3.54 (26) MeV, the pseudo scalar decay constant in the chiral limit f 0 = 122(1) MeV, the scalar condensate [Σ MS (µ = 2 GeV)] 1/3 = 270(7) MeV and f π /f 0 = 1.0755(94). The errors are statistical and systematical errors summed in quadrature. In case of asymmetric errors we use conservatively the maximum for the sum. The detailed error budget can be found in table 1 below.
A special emphasis of this paper is the investigation of systematic errors on the results for the low energy constants from the above mentioned chiral fits. To estimate the systematic errors, we perform different kind of fits where we take lattice spacing, finite size and next-to-next-leading order (NNLO) corrections into account. Finally, we address the question of isospin violations by comparing the neutral to charged pseudo scalar masses as well as discussing the isospin splittings for other physical quantities.
The paper is organised as follows: after describing the lattice set-up in the remainder of this section we shall present the main results of this paper in section 2. The following section 3 introduces the simulation set-up as well as the main ingredients entering in the continuum-limit scaling analysis of light meson observables. In section 4 we present the results of the combined chiral perturbation theory fits and section 5 gives a summary of our results. The details of the analyses and data tables are collected in the appendices.

Lattice Action
The twisted mass action, the tuning procedure to maximal twist and the analysis techniques as employed by our European Twisted Mass Collaboration (ETMC) have been discussed extensively in refs. [6,9]. We therefore only briefly recapitulate the essential ingredients of our set-up.
In the gauge sector we employ the tree-level Symanzik improved gauge action (tlSym) [29], viz.
with the bare inverse gauge coupling β = 6/g 2 0 , b 1 = −1/12 and b 0 = 1 − 8b 1 . The fermionic action for two flavours of twisted, mass degenerate quarks in the so called twisted basis [3,4] reads where m 0 is the untwisted bare quark mass, µ q is the bare twisted quark mass, τ 3 is the third Pauli matrix acting in flavour space and is the mass-less Wilson-Dirac operator. ∇ µ and ∇ * µ are the forward and backward gauge covariant difference operators, respectively. Twisted mass fermions are said to be at maximal twist if the bare untwisted quark mass m 0 is tuned to its critical value m crit , the situation we shall be interested in. For convenience we define the hopping parameter κ = 1/(8 + 2am 0 ). Note that we shall use the twisted basis throughout this paper.
Maximally twisted mass fermions provide important advantages over Wilson's originally proposed formulation of lattice QCD: the spectrum of Q † Q with Q = γ 5 (D[U] + m 0 + iµ q γ 5 ) is bounded from below, which was the original reason to consider twisted mass fermions [3]. At maximal twist, the twisted quark mass µ q is related directly to the physical quark mass and renormalises multiplicatively only. Many mixings under renormalisation are expected to be simplified [4,5]. And -most importantly -as was first shown in ref. [4], physical observables are automatically O(a) improved without the need to determine any operator-specific improvement coefficients. Another feature of maximally twisted mass fermions is that the pseudo scalar decay constant f PS does not need any renormalisation which allows for a very precise determination of this quantity.
The main drawback of maximally twisted mass fermions is that both parity and flavour symmetry are broken explicitly at non-zero values of the lattice spacing. However, it turns out that this is presumably only relevant for the mass of the neutral pseudo scalar meson (and kinematically related quantities). We shall discuss this issue in more detail later on. Note that in the following we use the notation m PS ≡ m ± PS for the charged pion mass, while the neutral one will be denoted by m 0 PS .
Tuning to maximal twist is achieved with the general prescription [30,31] to choose a parity odd operator O and determine am crit such that O has vanishing expectation value at fixed physical situation for all lattice spacings, e.g. by fixing m PS in physical units. One appropriate quantity to tune to zero is the PCAC quark mass, defined as where A a µ and P a are the axial vector current and the pseudo scalar density in the twisted basis, respectively, Once rotated to the physical basis the numerator of the r.h.s of eq. (2) represents the vacuum expectation value of a parity odd operator. Maximal twist is then achieved by demanding that the PCAC quark mass of eq. (2) vanishes. A discussion about the precise conditions we employ in our simulations is given in ref. [9] and summarised in section 3.

Main Results
In this section we present the main results of this paper. All the remaining sections are devoted to the discussion of the details leading to these results and to the estimation of systematic uncertainties.  extrapolating r 0 /a to the chiral limit. While our data cannot discriminate between these different ansätze, it turns out that all of our final physical results are independent of the exact way r 0 /a is chirally extrapolated. This is due to the fact that the physical scale is eventually set by the pion decay constant f π . The values of r χ 0 /a we used for the continuum-limit scaling analysis were obtained by chirally extrapolating r 0 /a linearly in (aµ q ) 2 at every value of the lattice spacing separately, as explained in section 3. Later on, in the context of combined chiral fits, we also allow r 0 /a to have an additional aµ q dependence, as discussed in section 4.
Leading lattice artefacts in our data are expected to be of order a 2 , as discussed previously. This can be checked by extrapolating a physical quantity at fixed physical situation to the continuum limit. We show two such examples in figure 1. In the left panel we show r χ 0 f PS as a function of (a/r χ 0 ) 2 at fixed value of r χ 0 m PS . In order to match the values of r χ 0 m PS at each value of r χ 0 /a and to fix the volume to r χ 0 · L = 5 we had to perform short inter-or extra-polations. The straight lines are linear fits in (a/r χ 0 ) 2 to the corresponding data, with the data at the largest value of the lattice spacing not being included in the fit. It is clearly visible that the lattice artefacts appear to scale linearly in a 2 and that their overall size is small.
In the right panel we show the scaling of r χ 0 m PS as a function of (a/r χ 0 ) 2 at fixed values of the renormalised quark mass r χ 0 µ R , again at fixed, but finite volume. We conclude that also the charged pseudo scalar mass has only small lattice artefacts. This observation will become particularly important when we shall discuss later the mass difference between charged and neutral pseudo scalar mesons. Note that we do not yet have the value of Z P for the smallest value of the lattice spacing which is hence left out in fig. 1 Table 1 Summary of fit results, determined from the weighted distribution as explained in the text. The first error is of statistical origin while the second, the asymmetric one, accounts for the systematic uncertainties. B 0 , Σ and m u,d are renormalised in the MS scheme at the renormalisation scale µ = 2 GeV, as the values of Z P are in the MS scheme at scale 2 GeV. Note that for the results listed here only data at β = 3.9 and β = 4.05 have been used. The scale is set by f π = 130.7 MeV as done in ref. [6]. For a comparison to other recent lattice results we refer the reader to ref. [2].
The dependence of m PS and f PS on the renormalised quark mass and volume can be described by chiral perturbation theory (χPT) [33][34][35]. The residual lattice artefacts of order a 2 can also be included in the analysis. The corresponding formulae can be found in equation (8) below. We fit these formulae to our data in order to extract the parameters of the N f = 2 chiral Lagrangian, i.e. the low energy constants and some derived quantities. Moreover, we can use these fits to calibrate our lattices by determining the value of the renormalised quark mass r χ 0 µ R where the ratio m PS /f PS assumes its physical value (i.e. m π /f π ) and set f PS = f π = 130.7 MeV there, 2 as done in ref. [6]. Hence, f π is used in this paper to set the scale.
The results of these fits can be found in table 1. We give statistical and systematic errors separately, the systematic one being asymmetrical. The results are obtained by performing O(80) fits, which differ in fit-range, finite size correction formulae and in the order of χPT. The final result is obtained as the median of the corresponding weighted distribution over all fits. The statistical error is determined using the bootstrap method with 1000 samples. The systematic uncertainty is estimated from the 68% confidence interval of the weighted distribution. For details we refer to the remaining sections of the paper, in particular to appendix A. For a comparison to other lattice results we refer the reader to ref. [2].
The results quoted in table 1 are obtained from the two intermediate values of the lattice spacing only, since at the smallest available lattice spacing Z P has not yet been computed and the largest value of a is not considered, since there our condition for tuning to maximal twist might not have been fulfilled accurately enough, which may affect in particular the observables f PS and m PS (at the smallest quark mass values) considered in this paper. However, including data at the smallest (by fitting Z µ = Z −1 P ) and the largest value of the lattice spacing gives results which are compatible within errors with the numbers shown in table 1, as can be seen from a direct comparison in table 7.
3 Ensemble details, r 0 , a 2 and finite size effects Before discussing the chiral perturbation theory fits to m PS and f PS in section 4 in detail we shall address in this section a number of issues which are needed as preparatory steps. These are the determination of the hadronic scale r 0 /a, the tuning to maximal twist, the size of the O(a 2 ) corrections (including isospin breaking effects) and finite size effects.

Ensemble Details
In table 2 we summarise the various N f = 2 ensembles generated by the ETM collaboration. We remark here that all the configurations produced by ETMC are available on the international lattice data grid (ILDG) (see ref. [37] and references therein), and the configurations are made public after publication of this paper. We have simulations at four different values of the inverse gauge coupling β = 3.8, β = 3.9, β = 4.05 and β = 4.2. The corresponding values of the lattice spacing, obtained from chiral fits as detailed later, are a ≈ 0.1 fm, a ≈ 0.079 fm, a ≈ 0.063 fm and a ≈ 0.051 fm, respectively. For each value of β we have several values of the bare twisted mass parameter aµ q , chosen such that the ensembles cover a range of pseudo scalar masses between 280 and 650 MeV.
The physical box length L of most of the simulations at β = 3.9, β = 4.05 are roughly equal and around L ≈ 2 fm, while the volume at β = 3.8 is slightly larger. At β = 4.2, the lattice size ranges from L ≈ 1.7 fm to 2.5 fm. For all the values of β, we demand the condition m PS L 3. Note that we have carried out several simulations at different physical volumes for otherwise fixed parameters in order to study finite size effects (FSE).  Table 2 Summary of ensembles generated by ETMC. We give the lattice volume L 3 × T and the values of the inverse coupling β, the twisted mass parameter aµ q , the critical hopping parameter κ crit as determined at µ q,min (i.e. the value of aµ q appearing in the same row as κ crit ) and the trajectory length τ . The values of the lattice spacing that correspond to the four values of β are a ≈ 0.1 fm (β = 3.8), a ≈ 0.079 fm (β = 3.9), a ≈ 0.063 fm (β = 4.05) and a ≈ 0.051 fm (β = 4.2). In addition we provide values for the integrated autocorrelation time of two typical quantities, the plaquette P and the pseudo scalar mass am PS , in units of τ = 0.5. We refer to ref. [9] for details on the determination of the autocorrelation time.
The simulation algorithm used to generate these ensembles is a Hybrid Monte Carlo algorithm with multiple time scales and mass preconditioning. It is described in detail in ref. [38] and one implementation described in ref. [39] is freely available. Its performance for the current set-up is discussed in ref. [26]. In table 2 we provide the values of the actual trajectory length τ used in the simulations. For each value of β and aµ q we have produced around 5000 equilibrated trajectories in units of τ = 0.5. In all cases we allowed for at least 1500 trajectories for equilibration (again in units of τ = 0.5). Some of the ensembles have been produced with a trajectory length τ = 1 following the suggestions of ref. [40].

Hadronic scale parameter r 0
In order to be able to compare results at different values of the lattice spacing it is convenient to use the hadronic scale r 0 [32]. It is defined via the force between static quarks and can be measured to high accuracy in lattice QCD simulations. For details on our procedure to determine r 0 /a we refer to ref. [9]. Note that in this work we use f π to set the scale and employ r 0 only for a continuum limit scaling analysis. In other words r χ 0 is used only as an intermediate reference quantity, which is finally eliminated in favour of f π .
In fig. 2 we show r 0 /a as a function of (aµ q ) 2 for the examples of β = 3.9 and β = 4.05. The mass dependence appears to be rather weak and a linear fit in µ 2 q describes the data well as has been also discussed in refs. [6,9]. However, also a linear dependence on aµ q cannot be excluded -due to possible effects of spontaneously broken chiral symmetry in r 0 . This dependence describes in the case of β = 3.90 the data even better than the quadratic ansatz. We shall therefore conservatively include both terms in our final analysis, i.e. the chiral fits discussed in section 4.
For tuning to maximal twist as well as for the continuum-limit scaling analyses presented in this section -in particular in figures 1, 3 and 4 -we shall use the value of r χ 0 /a from a linear extrapolation in (aµ q ) 2 to the chiral limit. A linear extrapolation in (aµ q ) 2 yields r χ 0 /a = 4.47(6) at β = 3.8, r χ 0 /a = 5.25(2) at β = 3.9, r χ 0 /a = 6.61(2) at β = 4.05 and r χ 0 /a = 8.33(5) at β = 4.2. The statistical accuracy for r χ 0 /a is about 0.5%. It is important to notice that the ratio [r χ 0 /a(β = 3.9)]/[r χ 0 /a(β = 4.05)] is within the errors independent of the extrapolation procedure. This leads us to expect rather little influence of the r 0 /a extrapolation strategy on physical quantities.

Tuning to Maximal Twist
Let us briefly discuss our strategy to tune to maximal twist. It corresponds to tune am 0 to a critical value am crit such that the PCAC quark mass defined in eq. (2) vanishes. In particular, we determine the value of am crit at the lowest available value of aµ q,min ≪ aΛ QCD at each β-value. Furthermore, we demand that the ratio |Z A am PCAC /aµ q,min | 0.1 and also that its error ∆(Z A am PCAC /aµ q,min ) 0.1. These criteria and their justification for tuning to maximal twist have been discussed in ref. [9]. An additional condition is that the value of aµ q,min corresponds to a fixed physical situation at all values of the lattice spacing. We have chosen to use m PS ≃ 300MeV to define the value of aµ q,min . Considering β = 3.9, β = 4.05 and β = 4.2, it was possible to perform this tuning task with two or three tuning runs for each lattice spac- ing. Obeying these conditions leads to very small O(a 2 ) effects in the physical observables considered in this paper, with the only exception of the neutral pseudo scalar mass, which will be discussed below.
A special case are the simulations at β = 3.8. As can be seen in table 2 at the smallest value of aµ q = 0.006 we encounter very large autocorrelation times for the plaquette. Since the autocorrelation time for the PCAC quark mass is of the same order of magnitude as the one of the plaquette, at this value of aµ q it is questionable whether am PCAC was reliably determined since the effective statistics is small. Hence, it is unclear whether a tuning to maximal twist fulfilling our conditions has been achieved at β = 3.8. A clarification of this point would have required substantially larger statistics. As a consequence, we will leave out the data at β = 3.8 in our main analysis. The β = 3.8 data will only be discussed in the context of systematic uncertainties.
The measurements of m PCAC from all the ensembles are collected in Appendix C. In figure 3 we show the properly renormalised ratio of the PCAC quark mass m PCAC over the twisted mass µ q (the renormalisation factor of Z P cancels in this ratio) against the renormalised twisted mass, µ R = µ q /Z P , in units of r χ 0 , for the four values of β. The renormalisation factors Z P , taken here at a scale of 2 GeV, were determined using the RI'-MOM scheme [41,42]. Note that for β = 4.2 the renormalisation factors are not available yet. However, we have estimated their values by fitting Z P and by estimating the value of Z A = 0.75 from the β-dependence of its known values at smaller β. Since these renormalisation constants enter as O(1) multiplicative pre-factors, a small change in their values will not significantly change the data in figure 3.
Renormalised ratio of the PCAC quark mass over the twisted mass against the renormalised twisted mass µ R = µ q /Z P at the four values of β. The statistical uncertainties on Z P and Z A are not included. The data at β = 4.2 have been included by estimating the renormalisation constants as described in the text. At β = 3.8, the data for the lightest quark mass has not been included for the reasons explained in the text. The band indicates our condition for tuning to maximal twist, which is clearly achieved to a good precision. The arrow indicates the value of r χ 0 µ R where we tuned the PCAC mass to zero.
Therefore our main conclusion about the quality of the tuning will be unaltered if the actual numbers for the renormalisation constants should come out slightly differently in the final analysis.
As can be seen from fig. 3 our condition for tuning to maximal twist is satisfied with a good statistical precision. The fact that at the three largest values of β the renormalised quark mass computed at the tuning twisted mass parameter aµ q,min basically agree (this is illustrated by the arrow in fig. 3) shows in addition that the physical situation has indeed been fixed. At β = 3.8, due to the previously mentioned problem of the long autocorrelation time at the lowest quark mass, we include in fig. 3 only the ensembles for which a reliable estimate of the error on m PCAC was possible. At all β-values for values of aµ q > aµ q,min we observe in am PCAC (small) deviations from zero. This O(a) cut-off effect will modify only the O(a 2 ) lattice artefacts of physical observables.

O(a 2 ) effects
In this section, we shall discuss the lattice spacing effects on the decay constant f PS and the mass m PS of the (charged) pseudo scalar meson. Starting with f PS , we will use the data at β = 3.9, β = 4.05 and β = 4.2. In order to compare the data at these values of β, we have fixed three reference values of the pseudo scalar mass to r χ 0 m PS = 0.614, r χ 0 m PS = 0.90 and r χ 0 m PS = 1.10. The corresponding values for af PS were then obtained by small interpolations, which we have chosen to be linear in m 2 PS . In addition, we also corrected for the very small differences in the physical volume for the β values used such that all data were scaled to the same physical volume r χ 0 · L = 5.0. This has been achieved by applying the relevant formulae from χPT as detailed below.
As has already been discussed in section 2, in fig. 1(a) (cf. page 6) we show the pseudo scalar decay constant, brought to common reference points r χ 0 m PS and volumes, as a function of (a/r χ 0 ) 2 . In our earlier work, where we had only two values of the lattice spacing [26][27][28], we detected only very small lattice spacing artefacts when comparing β = 3.9 and β = 4.05. This observation is confirmed: in fig. 1(a) we now have added, for two values of aµ q , a smaller lattice spacing of a = 0.051 fm and we can clearly observe that the scaling is as expected linear in a 2 with only a small slope such that no large cutoff effects are visible. Note that also for the nucleon mass the scaling violations are very small (see refs. [10,43]).
We also show r χ 0 f PS for β = 3.8 in fig. 1(a) at least for those points where we are able to interpolate to the reference values of r χ 0 m PS . Although the linear continuum extrapolation is performed from the results at β = 3.9, β = 4.05 and β = 4.2, as can be seen, the data at β = 3.8 are consistent with this linear behaviour. This indicates that the lattice artefacts are indeed small.
In fig. 1(b) we show similarly the continuum extrapolation of (r χ 0 m PS ) 2 for three fixed values of the renormalised quark mass r χ 0 µ R = 0.045, 0.09, 0.13. In this figure we cannot include data at the smallest lattice spacing, due to the missing Z P -factor. Instead we include data from β = 3.80 in the continuum extrapolation to show that the scaling looks very flat. In conclusion also for (r χ 0 m PS ) 2 residual scaling violations are small.

Isospin Breaking Effects
A peculiar O(a 2 ) effect in the twisted mass formulation of lattice QCD is the breaking of isospin symmetry at any non-zero value of the lattice spacing. Investigations of these effects deserve therefore special attention. In the quenched approximation the difference (r χ 0 ) 2 ((m ± PS ) 2 − (m 0 PS ) 2 ) between the neutral m 0 PS and the charged m ± PS pseudo scalar masses, although fully compatible with the expected O(a 2 ) behaviour, turned out to be significant [23].
The determination of m 0 PS requires the computation of disconnected diagrams rendering its determination difficult. Nevertheless, following the techniques described in ref. [9], we were able to compute m 0 PS (and other neutral quan-  Table 3 Comparison of some selected quantities for which an isospin splitting can occur for twisted mass fermions. R O denotes the relative size of the splitting, see also ref. [24,25]. tities) with a reasonable precision. The results for m 0 PS and m 0 V can be found in table C.10 in the appendix.
We show in fig. 4 the mass splitting between the charged and neutral pseudo scalar mesons, i.e. ( as a function of (a/r χ 0 ) 2 . As can be seen, also for N f = 2 dynamical quark flavours the isospin breaking effects are large, although smaller than observed in the quenched approximation and of opposite sign, cf. ref. [23], meaning that in the case of dynamical quarks considered here the neutral pseudo scalar meson turns out to be lighter than the charged one. This observation is consistent with the observed first order phase transition [44,45] and with the corresponding scenario of ref. [46].
One very important observation is that the large cutoff effect in the mass difference between the neutral and charged pseudo scalar masses is dominated by discretisation effects in the neutral pseudo scalar mass. This is demonstrated in fig. 1(b) which shows the charged pseudo scalar mass as a function of (a/r χ 0 ) 2 for fixed values of r χ 0 µ R . As can be seen, the lattice spacing effects are very small, as discussed above.
Another observation is that, in the unitary sector, so far all other investigated differences between corresponding quantities affected by isospin violations are compatible with zero. In table 3  The ETM collaboration has been investigating the question of the large cutoff effects in the neutral pseudo-scalar mass [24,25,47] also theoretically. An  (430) MeV. The open circle is a larger physical volume. The lines are only to guide the eye and some points are slightly horizontally displaced for better visibility. Note that r χ 0 m ± PS was not held fixed for this plot, however, due to the large uncertainties on m 0 PS the picture should not significantly depend on this approximation.
analysisà la Symanzik of the charged and the neutral pseudo scalar meson masses leads to the formulae which show that the difference (m 0 PS ) 2 − (m ± PS ) 2 is given by the term proportional to ζ π (for a discussion in the framework of twisted mass Wilson χPT see ref. [48,49]). Here ζ π ≡ π 0 |L 6 |π 0 and L 6 is the dimension six term in the Symanzik effective Lagrangian.
The main result of the analysis of ref. [24,25,47] is that ζ π is a large number which in the vacuum saturation approximation can be estimated to be proportional to |Ĝ π | 2 , whereĜ π = 0|P 3 |π 0 . The latter matrix element is numerically large: one finds |Ĝ π | 2 /Λ 4 QCD around 20 − 25 [24,25,47]. This result provides a physical explanation for the large O(a 2 ) effect in the neutral pseudo scalar mass. Moreover, since it can be shown that ζ π appears only in the neutral pseudo scalar mass (and kinematically related quantities), one also finds a possible explanation of why all other splittings determined so far turn out to be small. We remark that in principle also a double insertion of the operator L 5 can contribute to the O(a 2 ) lattice artefacts. However, it turns out that these contributions are much smaller than those of L 6 discussed here [25].

Finite Size Effects in f PS and m PS
In the following, we shall focus on the treatment of finite size effects in f PS and m PS with the aim of assessing the applicability of an effective description to be used, in particular, to perform the very small volume corrections required to bring the simulation points to a common finite volume, as needed for the continuum-limit scaling analysis. Later on, in the context of the chiral fits described in section 4, we will apply the finite size corrections formulae to bring our continuum data to infinite volume and to check for systematic effects.
Our lattice data have been obtained with values of m PS L 3. It is generally believed that in such a situation the finite size effects are small and appear only at the per cent level. However, our simulation data are on a level of precision that such effects are certainly detectable and hence finite size effects can in general not be neglected. One possible way to control finite size effects is to use chiral perturbation theory (χPT) which provides an effective description for the dependence of physical quantities on the box size [50]. A completely independent analytical description of finite size effects (FSE) of particle masses is provided by the so called Lüscher formula developed in ref. [51].
Being interested in m PS and f PS -where the applicability of χPT should be best justified -we then have two effective descriptions of FSE at hand: (i) χPT to next-to-leading order (NLO) by Gasser and Leutwyler [50], which we shall denote with GL for short, and (ii) the extended asymptotic Lüscher formula of Colangelo, Dürr and Häfeli [52,53], which conveniently combines the Lüscher formula approach with χPT and which we will refer to as CDH in the following. Note that in pure χPT the finite size corrections to m PS are also known to 2-loop order [54]. However, it turns out that the difference between the CDH formula and 2-loop χPT is negligible [54].
Below, we compare the numerical data obtained by ETMC to the two different analytical descriptions of the FSE in m PS and f PS , which is possible since we have precise values for these quantities for a number of ensembles such as C 2 , C 5 and C 6 or B 1 and B 6 where all parameters but the lattice size are fixed. We remark that for all these ensembles m PS L ≥ 3 holds, such that the formulae from refs. [50,[52][53][54] should be applicable. This was also checked and discussed in ref. [26].
Starting with the GL ansatz, the finite size corrections for m PS and f PS read where g 1 is a known function [50] and the finite size corrections K GL m,f depend apart from L and m PS only on the unknown leading order low energy constant f 0 representing the pseudo scalar decay constant in the chiral limit (note that our normalisation is such that the physical value of the pseudo scalar decay constant is f π = 130.7 MeV).
The CDH formulae have a very similar form where the I (i) m,f can be written in terms of basic integrals as discussed in ref. [53], where also m(n) is tabulated. The authors of ref. [53] provide simplified formulae for the I (i) m,f , which we shall use in the following. Note that I (6) m is known while I The drawback of the CDH formulae compared to GL is that additional low energy parameters are needed as an input 3 , namely: Λ 1 , Λ 2 , Λ 3 and Λ 4 . However, there are estimates for all those parameters given in ref. [55] and compiled in the forml i = log(Λ 2 i /m 2 π ) in table 4 for convenience, which we shall rely on during the analysis presented in this section. In order to convert the estimates in table 4 into lattice units we conventionally used in this analysis r χ 0 = 0.45 fm. Note that the choice of scale here affects only the sub-leading contributions to the finite size corrections and has a negligible effect, as was checked also by direct inspection. We emphasise again that we are at this stage only interested in an effective description of FSE with the goal to understand whether the correction formulae are applicable in the regime of our simulation parameters.
While the expansion in ref. [55] was performed in terms of the squared pion mass, in ref. [56] the same formulae have been obtained expanding in the quark mass directly. We shall refer to the latter as CDH m and include it in the following tests.
We can now apply the different formulae to correct m PS and f PS for finite size effects for ensembles where we have different volumes available. The result il i 1 −0.4 ± 0.6 2 +4.3 ± 0.1 3 +2.9 ± 2.4 4 +4.4 ± 0.2 Table 4 Values forl i = log(Λ 2 i /m 2 π ) from table 2 of ref. [53] used here to analyse FSE. (7) 0.1354 (7) 0.1348 (7) 0.1350 (7) 0.1346 (7   can be found in tables 5 and 6. In the tables we show the corrected values at infinite volume from the different ways of treating finite size effects as indicated in the caption. To this end, we take the measured data at a given linear extent L and provide the corresponding correction to infinite volume. For larger volumes these extrapolations agree better and finally converge within the errors. By checking the stability of the values extrapolated to infinite volume from the different analytical formulae we can control the applicability of these corrections for a given lattice size.
All finite size correction formulae provide an appropriate framework to describe the observed finite volume effects. However, for ensembles B 1 and C 6 , corresponding to the smallest lattice extents at these two β values, GL seems to underestimate FSE in m PS , an observation which was also made in refs. [57,58]. For those cases the resummed Lüscher formula provides a better description.
In the chiral fits presented in section 4, however, we observe significantly larger values of χ 2 /d.o.f when GL is used, and hence we use CDH and CDH m in these fits.
In conclusion, these results make us confident that our simulations have eventually reached a regime of pseudo scalar masses and lattice volumes where χPT formulae can be used to estimate FSE for m PS and f PS when ensembles with m PS L 3 are used. In particular, this allows us to control the finite size effects for these mesonic quantities from our simulations. But it is clear that in particular the CDH formula is affected by large uncertainties, mainly stemming from the poorly known low energy constants, which are needed as input.
Changing their values in the range suggested in ref. [53], however, changes the estimated finite size effects maximally at the order of about 20% (of the corrections themselves) [9].
Note that in the above discussion we have only used continuum formulae to describe the FSE of our lattice data. However, we believe that this can be justified by the smallness of the lattice artefacts discussed above. In addition, we restate that at this stage we are only interested in an effective but good description of the FSE since our goal was to bring the simulation data from different lattice spacings to a common finite volume as needed to perform the scaling analysis. In all the fits performed later on, the finite size corrections are applied to the continuum results.

Combined Fits
The main goal of this section is to confront the data for f PS and m PS to χPT and eventually determine low energy constants (LEC's) of the chiral effective Lagrangian.
We shall proceed by presenting results from fits including a combined continuum, infinite-volume and chiral extrapolation of m PS and f PS for two β-values, β = 3.9 and β = 4.05. Our analysis will extend the results given in refs. [26,27] by incorporating data for the renormalisation constant Z P and the Sommer parameter r 0 /a in the fit. Details on the computation of Z P and r 0 /a can be found in refs. [9,42]. We shall include the correlation matrix of f PS and m PS into the fit (the correlations with r 0 /a and Z P are negligible) and also provide a detailed account of systematic uncertainties.

General Fit Formulae
The formulae we are going to use to describe the chiral and continuum behaviour of f PS , m PS in infinite volume read: where T NNLO m,f denote the continuum NNLO χPT terms and we have defined We use a normalisation such that f 0 = √ 2F 0 , i.e. f π = 130.7 MeV. For the finite size corrections we shall use the CDH and CDH m formulae as discussed in sub-section 3.5.
The mass and decay constant of the charged pion have been studied up to NLO [48,49,59] in the context of twisted mass chiral perturbation theory (tmχPT). The regime of quark masses and lattice spacings at which we have performed the simulations is such that µ R aΛ 2 QCD . In the associated power counting, at maximal twist the NLO tmχPT expressions for the charged pion mass and decay constant preserve their continuum form. The inclusion of the terms proportional to D m PS ,f PS parametrising the lattice artifacts in eq. (8) represents an effective way of including sub-leading discretisation effects appearing at NNLO. The data shown in fig. 1 (page 6) nicely support the µ R independence of D f PS ,m PS , which we assume in our fit ansatz (8).
r χ 0 is the value of the Sommer parameter in the chiral limit which is determined for every β-value, within the global fit, according to the formula with β independent, i.e. continuum parameters C 1 and C 2 . In addition to the expected µ 2 R dependence, we also include the term linear in µ R in order to account for possible effects of spontaneously broken chiral symmetry. Possible lattice artifacts of O(a 2 ) cannot be distinguished from the lattice artifacts included for r χ 0 m PS and r χ 0 f PS in eq. (8), and are therefore effectively described by the parameters D m PS ,f PS , whereas lattice artifacts of order a 2 µ R and a 2 µ 2 R appear only at higher order.
In order to work with a massless renormalisation scheme, the values of the pseudo scalar renormalisation constant Z P are determined in the chiral limit. They were determined using the RI'-MOM scheme at scale µ ′ = 1/a and the values can be found in table C.5. For the present fit all renormalisation constants Z P need to be run from scale µ ′ = 1/a to a common scale, which we have conventionally chosen to be µ = 2 GeV. This running can be performed perturbatively, but the scale µ ′ is to be determined. This task can be performed for a given set of fit parameters by setting f π = 130.7 MeV where the ratio m PS /f PS assumes its physical value. Hence, we include for Z P into the global fit. ζ(µ, µ ′ ) describes the running of Z P from scale µ ′ = 1/a to the scale µ = 2 GeV computed in perturbation theory to four loops in α s in ref. [60]. We use the value Λ MS = 0.2567 GeV [61]. The perturbative factor R(µ = 2 GeV) = 0.8178 matches the RI'-MOM to the MS scheme and is known to NNNLO. It should be noted that the scale µ ′ enters only logarithmically into ζ and the small changes in µ ′ therefore do make only a very small difference in the running of Z P . However, we include eq. (13) into the global fit in order to correctly account for the statistical uncertainty of Z RI ′ P (µ ′ = 1/a).
At NLO, i.e. setting T NNLO m,f ≡ 0, and neglecting finite size corrections for the moment, there are the following free parameters to be fitted to the data for af PS , am PS , r 0 /a and Z RI ′ P (µ ′ = 1/a): where we indicate with the notation {...} β that there is one parameter for each β-value. When also the NNLO expressions for m PS and f PS are included there are four more parameters to be fitted to the data: As already mentioned, finite size effects are corrected for by using the asymptotic formulae from CDH and CDH m , which is consistently included in the fit. However, both additionally depend on the low energy constants Λ 1,2 .
It will turn out that we do not have sufficient data to determine all those parameters from the fit. For this reason we shall add priors for the NNLO low energy constants where needed. A description of the procedure to compute the χ 2 can be found in appendix B.
As mentioned above, in our analysis, we find that for a given value of aµ q and β only the data for am PS and af PS are correlated and we use a correlated χ 2 to account for these. Estimates for the correlation between af PS and am PS can be found in the tables C.6-C.9 in appendix C.

Analysis strategy
By in-and excluding certain data points, we build different data-sets. For a given data-set we perform the following four different fits: In the above list of different fits, setting D m PS ,f PS ≡ 0 corresponds to a constant continuum extrapolation of f PS and m PS . All these fits are applied to every data-set summarised in table A.1, and are repeated with both CDH and CDH m , respectively. While we include or exclude m PS and f PS data corresponding to the data-set, for r 0 /a we have always used only the largest volume available, since an effective description of finite volume effects is not available for this quantity.
These different fits include the systematic uncertainties of lattice spacing artefacts, finite size effects, higher order contributions in chiral perturbation theory and the extrapolation of r 0 /a to the chiral limit. Our strategy is to build the distribution of all possible fits using the available data-sets, weighting with the corresponding confidence level. This procedure is analogous to what has been done in ref. [62]. In table 1 we give our best estimates for the fit parameters following this procedure. The best estimates and also quantities derived from those are obtained as the median of the weighted distribution from a total of 76 different fits. All fits are repeated for 1000 bootstrap samples which are used to estimate the statistical errors, while the systematic uncertainties are estimated from the 68% confidence interval of the weighted distribution, as explained in appendix A. As an example we show in figure 5 the result for fit type B on data-set 1, which is also further discussed in the appendix.
We exclude fits from the analysis which we consider as being not reasonable, i.e. where the NNLO fit gives a worse description at larger masses than the NLO fit. This is a sign for not having sufficient data for such a fit. In table 1 there are a few quantities which show a strong asymmetry in the estimated systematics. This uncertainty is coming from the fact, that the estimate from NLO and NNLO fits for this quantities differs substantially. However, in the final result the NLO fits have more weight leading to the observed asymmetries.
In order to obtain further confidence in the uncertainties we quote in table 1 we can use the data at β = 3.8 and β = 4.20. In table 7 we compare fits including  these data with the results from table 1. It is clearly visible that the overall uncertainty we quote for all the quantities comfortably covers the differences among the three presented fit averages. This fact makes us confident that we indeed have a good estimate of the systematic uncertainties in our results. Note that Z P (β = 4.2) is a fit parameter for the fits including data at β = 4.2.
We think that this strategy takes into account all systematic effects that may affect our data and which are accessible to us. This allows us to provide values for the fit parameters and derived quantities with controlled statistical and systematic errors. In particular, since in our fit procedure, finite size corrections are applied in the continuum, effects of the neutral pseudo scalar mass on the finite size effects are absent. The details for our fits can be found in appendix A. Note that the main uncertainty to our results is stemming from our estimates of systematic errors. A systematic effect we are not able to estimate is the missing strange and, possibly, charm quark in the sea.
Eventually it is interesting to investigate the influence of the way r 0 /a is chirally extrapolated on our final results. In table 1, the systematic effect on r χ 0 /a is negligible because of the little variation in the treatment of r 0 /a between the different fits. Looking at table 7 one might be surprised that the values for C 1 and C 2 are both compatible with zero, even if in figure 5(a) a clear mass dependence is visible. The reason for this is that with both parameters the fit is overdetermined and the boostrap samples for C 1 and C 2 are strongly anticorrelated (−0.99). However, since we are in this analysis not primarily interested in the estimates for C 1 and C 2 , but more in the correct estimate of systematics for other quantities, we keep both parameters in the fits.
In order to address this systematic effect, we have repeated all the analysis using i.e. setting C 1 = 0, as compared to eq. (12). In table 8 we compare the results obtained with eq. (14) with the ones presented in table 1. The comparison leads to the conclusion that only the values of r χ 0 /a are significantly changed, being also statistically more precise when eq. (14) is used. Thus, the systematic uncertainty stemming from the r 0 /a extrapolation to the chiral limit is negligible for physical quantities.
In the tables 1, 7 and 8 we also quote a number of quantities that are derived from the basic fit parameters. The values for the low energy constantsl i ,   Table 7 We compare fit results of fits including data from β = 3.8 (first data column) and data from β = 4.20 (third data column) with the result where only data from β = 3.90 and β = 4.05 have been used (second data column). Note that in the case of including data from β = 4.20 the value of Z MS P is determined from the fit. All averages are weighted averages with the confidence levels of the individual fits. The errors are statistical and systematical, added in quadrature. B 0 , Σ and m u,d are renormalised in the MS scheme at the renormalisation scale µ = 2 GeV, as the values of Z P are in the MS scheme at this scale µ. The scale is set by f π = 130.7 MeV as done in ref. [6].
are determined using the physical value of the charged pion mass m π , which is the usual convention. The scalar condensate can be determined from the fit parameters B 0 and f 0 in the following way:  Table 8 We compare the main results of table 1, where the chiral extrapolation of r 0 was performed according to eq. (12), with results obtained with the purely quadratic extrapolation of r 0 /a to the chiral limit in eq. (14). For further details see table 1. and the scalar radius at next to leading order is given by

Summary and Outlook
In this paper we have demonstrated that Wilson twisted mass fermions, when tuned to maximal twist, show the expected linear scaling in a 2 with an almost negligible coefficient for both the pseudo scalar decay constant f PS and mass m PS (see fig. 1 on page 6). A very similar behaviour is seen for the nucleon mass [10,43]. It is important to remark that tuning to maximal twist has been achieved in our work by demanding that the renormalised ratio Z A m PCAC /µ min q of the PCAC mass m PCAC over the minimal twisted mass µ min q employed in our simulations satisfies |Z A m PCAC /µ min q | < 0.1. The same condition is used for the corresponding error of this ratio.
The precise data for f PS and m PS we have obtained at several values of the lattice spacing, quark masses and volumes allowed for a detailed analysis of systematic effects originating from having to perform a continuum, infinite volume and chiral limit. In particular, in this way it became possible to confront our data to analytical predictions from chiral perturbation theory and to extract a number of low energy constants and derived quantities as listed in table 1, page 7. The main physical results we obtain from this analysis are the light quark mass m MS u,d (µ = 2 GeV) = 3.54 (26) MeV, the pseudo scalar decay constant in the chiral limit f 0 = 122(1) MeV, the scalar condensate [Σ MS (µ = 2 GeV)] 1/3 = 270(7) MeV and f π /f 0 = 1.0755(94).
By performing O(80) different fits, taking into account or leaving out terms describing e.g. lattice artefacts or NNLO chiral perturbation theory and including or excluding data sets, we could achieve a reliable estimate of systematic errors. The final numbers for the fit parameters and their errors were then derived from the median, the bootstrap samples and the 68% confidence level of the (weighted) distribution of fit parameters originating from the different fits. Some quantities such asl 3 = 3.50 (31),l 4 = 4.66 (33) or the scalar condensate belong to the most precise determinations available today. For a comparison to other lattice results we refer the reader to ref. [2,63,64].
As a particular effect appearing in the twisted mass formulation of lattice QCD, we determined the size of isospin violation in various quantities. In this analysis we found, in accordance with theoretical expectations [24], that only the neutral pseudo scalar mass is significantly affected by the isospin breaking.
We conclude that Wilson twisted mass fermions at maximal twist provide a lattice QCD formulation allowing for precise computations of quantities in the light meson sector of QCD. Clearly, the next step is now to include dynamical strange and charm degrees of freedom in our simulations and work in this direction is in progress [65,66].

Acknowledgments
The computer time for this project was made available to us by the John von Neumann-Institute for Computing on the JUMP and Jugene systems in Jülich and apeNEXT system in Zeuthen, by UKQCD on the QCDOC machine at Edinburgh, by INFN on the apeNEXT systems in Rome, by BSC on Mare-Nostrum in Barcelona (www.bsc.es), by the Leibniz Computer centre in Munich on the Altix system and by the computer resources made available by CNRS on the BlueGene system at GENCI-IDRIS Grant 2009-052271 and CCIN2P3 in Lyon. We thank these computer centres and their staff for all technical advice and help. On QCDOC we have made use of Chroma [67] and BAGEL [68] software and we thank members of UKQCD for assistance. For the analysis we used among others the R language for statistical comput-ing [69].
This work has been supported in part by the DFG Sonderforschungsbereich/ Transregio SFB/TR9-03 and the EU Integrated Infrastructure Initiative Hadron Physics (I3HP) under contract RII3-CT-2004-506078. We also thank the DEISA Consortium (co-funded by the EU, FP6 project 508830), for support within the DEISA Extreme Computing Initiative (www.deisa.org). G.C.R. and R.F. thank MIUR (Italy) for partial financial support under the contract PRIN04. V.G. and D.P. thank MEC (Spain) for partial financial support under grant FPA2005-00711.

Appendices A Details of fits and discussion of systematics
In this appendix we discuss the details of our data analysis. For a more comprehensive discussion see ref. [36]. For every fit we compute besides the usual χ 2 value also the associated confidence levels (CL) (sometimes called goodness of the fit) CL(q = χ 2 , n = dof) = 1 − P (q/2, n/2) , where P is the incomplete Gamma function which corresponds to the cumulative χ 2 -distribution function for the given number of degrees of freedom F χ 2 (q, n). F χ 2 (q, n) itself represents (in the frequentists approach to statistics) the probability of finding a value of χ 2 ≤ q. A confidence level of CL(q) = x indicates then the fraction of times the computed value of χ 2 is larger than q, even if the fitted model is correct.
Since we follow a Bayesian approach by including prior knowledge for some parameters into our fits, the interpretation of CL is more a Bayesian credibility level. Hence, we are not going to determine the most probably correct model on the basis of the values of the CL's, but we shall rather incorporate all different fits into the final result. The CL's of each fit are then used as weights, specifying the impact of a given fit on the final result.
In more detail, for obtaining the final results we first sort out all physically unreasonable fits, as will be described below. Then we generate the CL-weighted distribution of a given fit parameter or of a derived quantity θ over all the retained fits. The expectation value of θ is estimated as the median of the weighted distribution. By performing this procedure for every bootstrap sample we determine the statistical error on our estimate for θ. An estimate for the systematic uncertainty is provided by 68.4% confidence interval of the weighted distribution. The estimate for the final error is obtained by adding the statistical and the systematic error in quadrature.
These systematic errors cover effects from lattice artefacts, from different orders in χPT, finite size effects and fitting range of the quark masses. We present separately the results including data generated at β = 3.8 and β = 4.2, because they are on a different level of accuracy due to insufficient precision in tuning to maximal twist at β = 3.8 and due to the missing estimate of Z RI P at β = 4.2. The results for these averages can be found in table 7, in-and excluding β = 3.8 and β = 4.2, compared to the result using only β = 3.9 and β = 4.05 data. The total residual error we quote for our final results covers the difference to the results including β = 3.8 and β = 4.2 data.
As mentioned before, we vary the data-sets in order to probe the influence of the fit range and the finite volume effects on the χPT fit. The data-sets we use are compiled in table A.1. We use both CDH and CDH m to correct for finite size effects to further check for the influence of those on our results.
The priors for r 0 Λ 1,2 and k M,F are necessary to obtain stable fits. The priors for r 0 Λ 1,2 are taken from table 4 and k M,F = 0 ± 10 is sufficient.

NNLO Fits
The NNLO χPT formulae provide in general a better χ 2 /d.o.f. than the corresponding NLO formulae. However, it seems that the larger masses with r χ 0 µ R 0.14 are necessary to stabilise these fits. If these masses are left out from the fits then the NNLO fit curves strongly at large values of r χ 0 µ R and actually undershoots the data at larger masses. This can be seen in figure A.1, where we show in the left panel the data for r 0 f PS as a function of r χ 0 µ R and the best fit function for fit C. Note that the data at large values of r χ 0 µ R are not included in the fit. In the right panel we show the LO, NLO and the complete NNLO function for comparison. Evidently, for larger masses, there is a significant difference between the NLO and NNLO results. It is in particular suspicious that the NNLO fit shows a stronger curvature than the NLO fit and that the NNLO fit deviates more from the data points at large masses -which are not included in the fit. On the contrary one would expect the NNLO formulae to describe the data points at larger masses better than the NLO formulae. On the other hand, when including the heavier masses, the NNLO fit is able to describe these data points. To improve the sensitivity of our lattice data to χPT at NNLO, additional data points would be needed.
We remark that in ref. [56] the inclusion of NNLO χPT was found to be necessary in order to obtain an adequate fit to the data on the pion radius. In ref. [56] it was argued that due to the ρ-meson dominance one should expect a slower rate of convergence of chiral perturbation theory. However, using the final analysis is that the estimated systematic uncertainty in e.g. f 0 ,l 3 andl 4 is asymmetric. If we analyse NNLO fits only, we would compared to table 1 rather obtain f 0 = 122.6(1.1),l 3 = 3.22 (66),l 4 = 4.33 (34). This uncertainty is included in our final error estimate of those quantities.

Lattice Artifacts
When performing fit B, which uses only NLO χPT and includes lattice artifacts, we find that the fitting parameter D m PS parametrising the lattice artifacts in m PS is compatible with zero within errors. In f PS lattice artifacts are more apparent, even if the error on D f PS is still large.
However, from the fits we clearly observe that the χ 2 values are significantly smaller when lattice artifacts are included, and the weighted distribution is hence dominated by the fits of type B and D. As an example for such a fit we show in figure 5 plots for fit B on data-set 1. We show the data together with the fitted curves for (r 0 /a)/(r χ 0 /a) as a function of (r χ 0 µ R ) 2 in (a), for r χ 0 f PS in (b), for (r χ 0 m PS ) 2 in (c) and for (r χ 0 m PS ) 2 /(r χ 0 µ R ) in (d), the last three as a function of r χ 0 µ R . Note that in figures 5(b), 5(c) and 5(d) we did not include the error of r χ 0 and Z P for the data points.
From figure 5 one can observe that the fit works very well, leading to χ 2 /dof = 19/17. In (a) one can see that the data for r 0 /a is also fully compatible with a linear dependence on (r χ 0 µ R ) 2 only, as was also discussed earlier and summarised in table 8.

Finite Size Corrections and Effects of Priors
We have used the CDH and CDH m formulae for the fits and include both into the weighted distribution of all fits. It is still interesting to compare two fits, using CDH m on the one hand and the parameter free GL formulae on the other hand, while all the rest is identical. We observe that in general the χ 2 -value using CDH m is about a factor of 2 smaller than the one using GL. Therefore, we think that having extra parameters (which are included in the fit with priors) is justified.
Another test for uncertainties stemming from FS corrections is to perform fits where only the largest available volume per µ q -value is used. Including these ensembles (which correspond to data-sets 6,7 and 8) tend in general to produce improved values of χ 2 , while the influence on other fit parameters is again small. Clearly, it would be desirable to have for all ensembles values of m PS L > 4, where all the FS correction formulae give very similar results and the uncertainty from finite volumes is hence much reduced. In this investigation we include this significant uncertainty in our final error estimate.
As discussed previously, we have to add priors forl 1,2 and k M,F in order to stabilise the fits. It is therefore interesting to discuss the influence of the choice for the priors on the fit results. First of all, the influence of the choice for k M,F and their errors has negligible effect on our fit results.
In case ofl 1,2 the influence depends on whether or not NNLO terms are included in the chiral expansion. While doubling the error estimates ofl 1,2 or changing their values within the error range has negligible effects in the case of a NLO chiral fit (fits A and B), many fit results change drastically in the case of NNLO chiral fits (fits C and D): we observe in the latter case that the best fit values e.g. forl 3 andl 4 differ substantially from the NLO fit result and also from the estimates given in table 4. At the same time the errors become much larger on those quantities, up to 50%.
Moreover, the χ 2 minimisation process becomes more difficult, the χ 2 functions shows many nearby local minima in parameter space. These facts confirm once more that f PS and m PS data are not sensitive to NNLO terms in the chiral expansion.
Chiral Extrapolation of Z RI ′ P The chiral extrapolation of Z RI ′ P is not performed in the fits, but the chiral value plus error estimates is input to the fit. In order to check whether this might lead the fits into a wrong direction, we have performed fits with doubled uncertainty on the values of Z RI ′ P . It turns out that this has negligible effect on our results.

B χ 2 Implementation
We are attempting to fit partially correlated data including priors for some of the parameters. With our method we closely follow the approach discussed in ref. [70].
The data for f PS , m PS , r 0 /a and Z RI ′ P at the various values of the lattice spacing and the different volumes is described by the following set of fit parameters r χ 0 f 0 , r χ 0 B 0 , r χ 0 Λ 1−4 , k M,F , C 1,2 , {r χ 0 /a} β , {Z MS P (2 GeV)} β , D m PS , D f PS .
With the notation {...} β we mean that there is one fit parameter for each lattice spacing involved in the fit. Those parameters can be used, employing eqs. (8,12) and (13), to compute predictions for am PS , af PS , r 0 /a and Z MS P . We shall denote these predictions by P m PS , P f PS , P r 0 , P Z P and they all depend on a sub-list of fit parameters and on the spatial extend L/a, as presented in section 4.1.
We adopt the following definition for the χ 2 -function at fixed value of β: (B.1) The matrix C is the properly normalised covariance matrix, see e.g. refs. [71,72]. The last two terms in eq. (B.1) are included in order to appropriately include the measurements for r 0 /a and Z RI ′ P . Values for the correlation among m PS and f PS are given in the tables C.6-C.9. The full χ 2 is obtained by summing χ 2 (β) over all β-values.
In case we need to add prior knowledge for a fit parameter θ in order to stabilise the fit we include additional terms of the form to the fit, assuming Gaussian error distribution for θ with standard deviation δ[θ]. Here the predicted values P θ are the fit parameters and the "data" with errors are the priors. The term eq. (B.2) follows from Bayesian statistics, see ref. [70] and references therein.
In detail the computation of χ 2 is performed in the following steps for fixed values of the fit parameters: (1) compute infinite volume and continuum predictions: compute the continuum predictions for m PS and f PS using the formulae (8) in infinite volume for given values of r χ 0 µ R . Determine r χ 0 /a for each β-value according to eq. (12) and use it to estimate the scales 1/a from f π .
(2) compute finite volume predictions in the continuum: now apply finite size corrections to scale the data to the finite volumes used in the lattice simulations. the finite volume and finite a predictions for m PS and f PS can now be compared to the data and the χ 2 can be computed as detailed above. The χ 2 includes also the corresponding terms for the chiral extrapolation of r 0 /a and for the estimate of Z MS P .
The errors are estimated using the bootstrap method. For m PS and f PS we use bootstrap samples as generated from the raw data. For r 0 /a and Z RI ′ P those are generated from a Gaussian distribution with mean and standard-deviation properly chosen. This is justified since we checked that the correlation of m PS and f PS to r 0 /a and Z RI ′ P is negligible. This procedure allows to estimate errors on primary quantities, such as fit parameters for instance, as well as on derived quantities, like e.g. the lattice spacing.