A new method for fitting the low-energy constants in chiral perturbation theory

A new set of the next-to-leading order (NLO) and the next-to-next-to-leading order (NNLO) low-energy constants $L_i^r$ and $C_i^r$ in chiral perturbation theory is obtained. These values are computed using the new experimental data with a new calculation method. This method combines the traditional global fit and Monte Carlo method together. The higher order contributions are estimated with this method. The theoretical values of observables provide a good convergence at each chiral dimension, except the NNLO values of the $\pi K$ scattering lengths $a_0^{3/2}$ and $a_0^{1/2}$. The fitted values for $L_i^r$ at NLO are closed to their results with the new method at NNLO, i.e. these $L_i^r$ are nearly order-independent in this method. The estimated values for $C_i^r$ are consistent with those in the other literature as far as possible. Their possible upper or/and lower boundaries are also given. The values of some linear combinations of $C_i^r$ are also given. They are more reliable. If one knows the more exact values of some $C_i^r$, some other $C_i^r$ can be obtained by these values.

the accuracy of experimental data too. So far, a lot of researches arise and some LECs have been fitted. L r 1 , L r 2 and L r 3 are obtained by fitting K 4 form factors and ππ scattering lengths [6]. L r i (i = 1, 2, 3,5,7,8) are obtained by fitting the quark-mass ratio m s /m, the decay constant ratio F K /F π and K 4 form factors [7]. About ten years later, Ref. [8] adds ππ scattering lengths (a 0 0 and a 2 0 ), πK scattering lengths (a 1/2 0 and a 1/3 0 ) and the threshold parameters of the scalar form factor ( r 2 π S and c π S ) in the fit. Recently, Ref. [9] adds two-flavour LECs in the fit. The last two references not only obtain the next-to-leading order (NLO, O(p 4 )) LECs L r i , but also estimate a part of LECs C r i in the next-to-next-to-leading order (NNLO, O(p 6 )). Nevertheless, their results only make use of the theoretical expansion to a finite order (NNLO). The higher-order contributions are ignored. Some other problems also arise: (i) Some NLO fitted values of L r i have much difference with those obtained at NNLO. For example, L r 1 = 0.53(06)× 10 −3 at NNLO fit, which is about half of its NLO fitted value L r 1 = 1.0(1) × 10 −3 . Nevertheless, L r 1 is a constant. Its value is independent on the fitting chiral dimension.
(ii) The higher-order effects are not taken into account. If the higher-order corrections are considered, some physical quantities may change largely, such as a 1/2 0 and a 3/2 0 . Without this consideration, one can not know which physical quantities are fitted well. Hence, the errors of LECs can not be estimated more precisely.
(iii) The fitted χ 2 /dof is approximately equal to 1.0/10 at NNLO fit. There seems meet an overfitting problem.
Appropriate large χ 2 can give a wider range of C r i . They may lead to a better convergence of observables.
(iv) The πK scattering lengths a (v) Some original data about C r i [21] are very closed to the final results. The differences are less than 10 −12 . We guess these LECs are dependent on the boundaries.
In this paper, we attempt to solve the first three problems and study why the other two problems arise. We assume that the physical quantities are convergence by chiral dimension; L r i are nearly unchanged between the NLO and the NNLO calculations; C r i are not far away from those in the other references. A more reasonable set of L r i is obtained, and the estimates of C r i are also given. With some numerical calculations, the possible origins of the last two problems are also explained. This paper is organized as follows: In Sec. II, some hypothesises are given and the following calculations are based on these hypothesises. In Sec. III, the experimental data is given. Section IV introduces a modified global fit method with the higher-order estimation. In Sec. V, some rough values of L r i are given and the convergences of some observables are also presented. In Sec. VI, a new method is introduced, which can compute more reasonable values of L r i and estimate the values of C r i . Section VII gives the results of L r i and C r i with this new method. A short summary is given in Sec. VIII

II. THE LOW-ENERGY CONSTANTS AND THEIR HYPOTHESES
In ChPT, without the contact terms, for the three-flavor case, there are 10 LECs L i at NLO and 90 LECs C i at NNLO; for the two-flavor case, 7 LECs l i exist at NLO and 52 LECs c i exist at NNLO. Their renormalized values L r i , C r i , l r i and c r i are defined in Refs. [2,3,9,[22][23][24]. Some scale-independentl i [3] are used frequently. This paper will determine the values of 8 L r i (i = 1, . . . , 8) and 38 C r i (the values of i can be found in Table V) at the renormalization scale µ = 0.77GeV. Fourl i (i = 1, 2, 3, 4) will be used in the estimation and none of c i will be used. The following notations are the same as those in Ref. [9].
However, due to the experimental condition, the experimental data is lack. Until now, only 17 observables have been used in the global fit [9]. Theoretically, it is impossible to obtain all LECs very accurately by these observables. Hence, the precisions of these two order LECs are required differently. For L r i , their values are required more precise and more reliable, because their number (8) is much less than 17. On the other hand, although 17 is less than 38, it does not seem too small to estimate C r i . The intervals of C r i can be limited to some reasonable ranges at the least. To achieve these goats, the following hypothesises are introduced to limit the feasible ranges of the LECs.
(i) Chiral expansion for most observables is all assumed to have a good convergence. It means that the contribution from the high order is less enough than the low-order one. This is a theoretical assumption in the effective theory. Now, observables are expanded by the momentum and the quark masses in ChPT. For most observables, the LO values give most contributions. The NLO and the NNLO values are smaller and smaller. The sum of the unknown higher-order contributions, which is also called truncation error, should be smaller than the NNLO values.
(ii) All L r i are assumed to be stable. In other words, the values of L r i obtained at both NLO and NNLO are almost unchanged. This is because that LECs are all constants, and they are independent on the different computational methods. According to Hypothesis i, the contributions at NNLO and the truncation error would be small enough, and these small contributions only lead to a small variation of L r i . However, this does not always work. In Ref. [9], some NNLO fitted L r i have large difference from the NLO fitted ones. For example, the NNLO fitted value of L r 1 = 0.53(04) × 10 −3 is half of its NLO fitted value 1.0(1) × 10 −3 . The deviation from the NLO value is about 5σ. The reason is that these values are only fitted at a given order and the truncation error is neglected. For some L r i , this effect is not obvious, but for the other ones, it may have a large impact on their values. This hypothesis will be used for constraining the ranges of L r i at NNLO calculation, which needs be closed to the NLO calculation.
(iii) All C r i should be consistent with those obtained from the other references, i.e. their values can not deviate too much from those in the references. C r i do not like L r i , their number is large and their values are rare in the literature. Appendix B presents all relative C r i we can find. Their values distribute over some wide ranges, and not all of them agree with each other. The possible values of C r i are limited to these wide ranges. However, according to the first two hypotheses, some intervals for L r i can be excluded by these wide ranges.

III. OBSERVABLES, INPUTS AND χ 2
This paper is based on Refs. [8,9], which adopt a global fit method to obtain L r i and use a random walk algorithm to estimate C r i . For the NLO fit, the following 12 observables are used. The mass ratio m s /m can be calculated according to pion and kaon masses or pion and eta masses [7,8,17,21], where m s is the strange quark mass and m = (m u + m d )/2 is the isospin doublet quark mass. The ratio of the kaon decay constant F K to the pion decay constant F π (F K /F π ) is also used in the fit [8,9,17,21], which eliminates the unknown constant F 0 . There exist two form factors F and G in K 4 decay, their values and slops at threshold (f s , g p , f s and g s ) [7] are also considered in the fit. The ππ scattering lengths a 0 0 and a 2 0 [21,25], the πK scattering lengths a 1/2 0 and a 3/2 0 , and the pion scalar radius r 2 π S in the form factor F π S (t) [26] are also included. With these 12 observables, eight L r i (i = 1, . . . , 8) will be fitted. The other five observables are added at the NNLO fit, they are the pion scalar curvature c π S [26] and four two-flavour LECsl i (i = 1, . . . , 4) [27]. These observables are related to 8 L r i and 38 C r i . The total number 46 is larger than the number of observables. We will use a different method to obtain them. The renormalization scale µ is chosen to be 0.77 GeV in this paper.
The values of the meson masses and the pion decay constant are The average kaon mass is which is used in the calculation for the pion and kaon decay constants and the pseudoscalar meson masses [28]. The values of m s /m and F K /F π are [29] For K 4 form factors F and G, their slope and value at threshold are [29] f s = 5.712 ± 0.032, f s = 0.868 ± 0.049, g p = 4.958 ± 0.085, g p = 0.508 ± 0.122.
The latest results for π π scattering lengths are given in Ref. [30], which are based on the analysis of K e4 data. Their values are For π K scattering lengths, Ref. [31] gives the most recent experimental value for S-wave isospin-odd πK scattering length a − 0 = |a separately. Hence, we still use the same data as those in Ref. [9,32] Since no update has been found, the scalar radius r 2 π S and the scaler curvature c π S of the pion scalar form factor are the same as those in Ref. [9]. Their values are based on the dispersion analysis [33,34], For two-flavour LECsl i (i = 1, . . . , 4), the values ofl 1 andl 2 are chosen [35], which are the same as those in Ref. [9]. Forl 3 andl 4 , Ref. [9] uses the average of Lattice results [36,37] and the continuum results [3,35]. At that time, the Lattice results in Ref. [37] are not included in FLAG average [36]. The most recent FLAG [38] provide the following averages The values in Eq. (9) have included the results in Ref. [37]. A new estimate according to Eq. (9) and Ref. [3,35] is These are all physical quantities used in our calculation. The objective function in the estimation, χ 2 , is the same as those in Refs. [7][8][9] where X i(exp) are the experimental values, X i(th) are the theoretical estimates and ∆X i are the experimental errors. Generally, χ 2 is as small as possible. This function is a criterion to judge whether the LECs are reasonable or not.

IV. METHOD I: A MODIFIED GLOBAL FIT FOR OBTAINING L r i
In this section, a modified global fit method is introduced, which contained the higher-order estimates in the fit. This method is only estimating the values of L r i . The calculating process is similar to that in Refs. [8,9]. Only the differences are explained.

A. Chiral expansions
In ChPT, physical quantities are calculated order by order, but some quantities described above are mixed by different orders. In order to pick out the exact contributions from different orders, they need to be expanded order by order.
The expansion for the ratio F K /F π to the NNLO is [9] Hereafter, the subscript 2, 4, 6 and 8 are represented the contribution at the LO, NLO, NNLO and NNNLO, respectively.
The quark-mass ratio m s /m can be calculated according to the LO pion and kaon masses or the LO pion and eta masses Their expansions are (15) r 2 π S and c π S are related to the differential of form factors F π S (t). Their expansions are The first terms in Eqs. (16) and (17) are both equal to zero, it is due to the scalar form factor at LO is independent on t.

B. The estimation at the higher order
In the previous fitting methods [7][8][9], the influences from the higher orders have not been taken into account. Although the contributions from the higher orders should be very small according to Hypothesis i, it is also worth to evaluate the influences from the higher orders according to Hypothesis ii. Higher-order contributions may have a big impact on some values of L r i . However, the contributions of the order beyond NNLO are absolutely unknown, and they need to be estimated in other ways. Refs. [39,40] provide a method for quantitative estimation of the truncation errors, which is based on Bayesian method. They assume that the expansion coefficients of observables in the effective field theory are in natural size, and their distributions are symmetric about the origin. The distribution of the truncation errors is also symmetric about the origin. Their confidence intervals can be obtained by some ways. This assumption leads to zero truncation error. In practice, the contributions from the higher orders may not be equal to zero. Some non-zero estimates need to be obtained. In this section, a method for estimating the higher-order contributions is introduced, which is similar to that in Ref. [41], but we do some simplifications for saving computation time.
In ChPT, physical quantities are calculated by chiral dimension. Each order provides a small factor = p/Λ χ . For example, a physical quantity X can be written as where the dimensionless parameter Q = 2 , c n are dimensionless coefficients and X ref is the natural size of X. We take X ref equal to the LO value of X.
In practice, strict calculations in the higher orders are very complex because of a lot of unknown LECs and loop diagrams. Hence, the expansion of X is truncated at a certain order and only the first few terms can be obtained. If X is truncated at order k, the theoretical prediction for observable X is where k = 1, 2 and 3 represent the LO, NLO and NNLO, respectively. The truncation error is Before fitting the LECs, a non-zero and valid value of ∆ k needs be estimated first. For convenience, we estimate X directly, but not ∆ k or X . According to Hypothesis i, the sequence {c n Q n } is naively assumed to be a geometric sequence {a 0 q n }. Whether this assumption is reasonable or not depends on the final fitting results. This will be mentioned later. In this case, one can get the approximation, and the truncated error ∆ k is In order to determine the parameters a 0 and q in the geometric sequence, we definite two cumulative sums sequences where the sequence {S k } can be regarded as a set of discrete data, and they can be calculated if a set of LECs related to X is known. The cumulative sum S * k of the geometric series is The parameters a 0 and q can be fitted by the least squares method. For the NLO fit, only the LO and the NLO contributions of F π S (t) can be calculated. Then only the NLO contributions of r 2 π S and c π S can be calculated according to Eqs. (16) and (17). The higher-order contributions can not be estimated, because the LO contribution is zero. Hence, we assume F π S (t) as a geometric series, then r 2 π S and c π S are calculated simultaneously. For K 4 form factor, f s and g p are also the derivative of the form factor F and G, respectively. Their calculations are the same as those for r 2 π S and c π S .

C. Convergence
In order to give more constraints on LECs, besides the observables F K /F π , m s /m| 1 and m s /m| 2 , some other physical observables, i.e. F π , F K , m π , m K and m η are also introduced separately. We find that not all of them have good convergence simultaneously. Sometimes, the two-flavour LECs (l r i ) exp (i = 2, 3) also have bad convergence. The NNLO contributions of these observables may be larger than the NLO ones. Hence, we add the following new constraints to χ 2 as Ref. [9], where the denominator ∆ = 0.1, and f χ (x) = 2x 4 /(1 + x 2 ). Ref. [9] only adopts the first equation.  Table I presents the NLO fitted L r i . The results in the second column (Fit 1) assume L r 4 ≡ 0 and the other L r i are left free in the fit. The results in the third column (Fit 2) are obtained by a free fit. To compare with our results, the fourth (fifth) column lists the results in Ref. [9], which are fitted at NLO (NNLO). Both Fit 1 and Fit 2 are very closed to the NNLO fit in Ref. [9], but some of them are very different from its NLO fit. It indicates that the geometric series can give a good estimate for higher-order contributions. The estimation in Sec. IV B is valid. For Fit 2, L r 4 is small enough and it is satisfied the large-N c limit. This is an assumption in Ref. [9]. Its averages three sets of results by assuming L 4 = 0, ±0.3 to obtain the NLO results. Moreover, 2L r 1 − L r 2 and L r 6 are also satisfied the large-N c limit. They are also better than those in Column 4. It means that the estimates from the higher order can not be ignored. They have a great influence on L r i (especially L r 1 , L r 3 , L r 4 and L r 6 ). The large-N c limit is satisfied automatically. Hence, we have a good reason to believe that contributions beyond the NNLO also have a great influence on C r i . When we calculate the NNLO LECs, the truncation errors need to be estimated too. Since Fit 2 has no assumption about L r 4 and its values are not very different from Fit 1, we use it as the NLO results in this section.
The second to the fourth column in Table II shows the LO, the NLO and the higher-order contributions of the observables with Fit 2 L r i in Table I. The fifth column is the theoretical estimates. In order to see the convergence of these quantities obviously, the percentage of each order is defined, where X th is the theoretical estimate and the subscript "order" represents LO, NLO and higher order (HO). These percentages are shown in the parentheses in the second to the fourth columns in Table II. It shows that all observables have a good convergence. The χ 2 i = 2.1 from π K scattering lengths (a are not very small. This is the reason why there are large deviations between the fourth column and the fifth column in Table I, such as L r 1 , L r 3 , L r 4 and L r 6 . The convergences of m s /m| 1 , m s /m| 2 and F K /F π are already presented in Table II. If ChPT has a good convergence, their numerators and denominators also need to be convergent separately. The NLO contributions of them are where the values in the parentheses in the first row are (m 2 α ) 4 /m 2 α (α = π, K, η). These observables also have satisfactory convergences. In the NNLO fit, the greatest difficulty is that 46 unknown LECs C r i are involved. At present, we only find two methods can obtain all of their values. The latest resutls are in Refs. [9,16]. This subsection adopts these two sets of values first. We will estimate them in Sec. VI.
There is a little different from the NLO fit. Sec. V A has mentioned that πK scattering lengths a with the method in Sec. IV B. For example, the LO contribution of a 3/2 0 m π is −0.0709 in Talbe II, but the experimental value is −0.0448. If the NLO contribution has a small positive value, the NNLO contribution needs a larger positive value. Nevertheless, the first and the third term (a 0 and a 0 q 2 in Eq. (21)) in a geometric sequence have the same sign. The NNLO contribution is negative too. Hence, we assume that a In this case, ∆ a 1/2 0 is small if q 1 is small. The estimation of ∆ a 3/2 0 is the same as ∆ a 1/2 0 . For two-flavor LECs, Ref. [27] gives the relation betweens l r i , L r i and C r i up to the NNLO. In the NNLO fit, l r i (i = 1, . . . , 4) are estimated first, which depend on the renormalized scale µ. Thenl i are calculated by the estimates of l r i . The NNLO fitted L r i are shown in Table III. Column 2 and 4 use the C r i in Refs. [9] and [16], respectively. Column 3 and 5 are the relative deviations of L r To compare with our results, the results in the last column are the NNLO fit in Ref. [9]. For Fit 3, χ 2 /dof= 14.7/9 seems a little large. The main problem is that some Pct L r i in Column 3 are larger than 20%, such as Pct L r 1 , Pct L r 6 and Pct L r 7 . We consider these deviations are a little large. The value less than 20% is acceptable. The results for Fit 4 are even worse. χ 2 /dof= 79.8/9 is very large and most of Pct L r i are larger than 20%. Especially, the values of L r 2 and L r 7 are very different from their NLO fitting results. It indicates that these two sets of C r i in the references can not fit the data well at NNLO. A new set of C r i needs to be found. It needs satisfy all hypotheses in Sec. II.

VI. METHOD II
This section gives a new method to obtain a better set of L r i and C r i simultaneously. A part of processes in this method is similar to those in Method I.
Refs. [8,9] estimates C r i first, with a random-walk method in some restricted spaces of C r i . Then they fit L r i with the values of C r i . Although this method attempts to restrict the fitting values of L r i , some values of L r i still deviate too much from their NLO fitting values (see Table I). For example, L r 1 = 0.53(06) × 10 −3 and L r 2 = 0.81(04) × 10 −3 at the NNLO fit, which are about half of their values at the NLO fit L r 1 = 1.0(1) × 10 −3 and L r 2 = 1.6(2) × 10 −3 . We attempt to determine C r i by scattering points randomly in an area, but we do not find a better set of C r i . In addition, we find that C r i may not satisfy the large-N c limit. If we choose the none-zero C r i in Ref. [16], but the zero C r i are replaced by those in Ref. [9], a much smaller χ 2 can be found. Hence, we do not limit C r i too much. Fit 2 in Table I is used as a reference. According to Hypothesis ii in Sec. II, Pct L r i is assumed less than 20%, in other words, L r i are limited to The total number of L r i and C r i is 8 + 38 = 46, which is much larger than the number of the observables 17. There are 29 redundant parameters, and they can not be obtained exactly simultaneously. Hence, we expect all L r i are as precise as possible, but for C r i , large errors are possible. We adopt the following steps to obtain all of them. (i) All L r i are generated randomly according to uniform distribution in the ranges in Eq. (32). Before any fit, none of C r i is known and none of the theoretical values for the observables is known too. If C r i are adjusted to a reasonable set of values, any set of L r i may give a good fit. One can not eliminate any one of them at the beginning. Hence, the weight of each set is assumed to be equal, so uniform distribution is chosen. Actually, it is not always true. Some sets of L r i could not find a valid set of C r i . The related χ 2 /dof may be large, because the other hypotheses in Sec. II may not be satisfied. Hence, Some sets in Eq. (32) may be picked up. In this step, 6.8 × 10 5 sets of L r i are generated. This number is large enough, in order to keep enough sets of L r i at last. (ii) To avoid the unnecessary calculation, before obtaining C r i , some sets of L r i can be removed first. Most sets of L r i can not give a good prediction for l r i , i.e. the NLO theoretical values (l r i ) 4 deviate too much from their experimental values (l r i ) exp . These sets of L r i can not satisfy Eq. (27) and they can be removed. We use the following constraints first, The above constraint seems very weak, but most set of L i are constrained. After this steep, only about 6.8 × 10 4 sets of L r i are left. (iii) For a given set of L r i , the number of redundant parameters is 38 − 17 = 21, which is still large. One can not obtain a unique set of C r i . The random-walk method [8,9] may give a reasonable set of C r i , but the efficiency is low and the obtained values may be not always equal. Hence, We do not determine C r i directly. There exist 17 observables. 17 combinations of C r i can be determined uniquely. Generally, not all of these combinations are linear. Appendix A gives a method to change them to linear ones called C i . These C i are considered as a whole and they are computed first.
In this step, only the constraint in Eq. (27) is used, because A i (i = 1, . . . , 5) (defined in Appendix A) can not be obtained separately. The NNLO values of m 2 α (α = π, K, η) and F α (α = π, K) can not be obtained too. They are going to be put back later.
The values of C r 14 , C r 15 and C r 17 can be directly obtained by solving the linear equations where P ij is a coefficient matrix, j is not continuous and its value is belonging to the numbers in the first column in Table VIII. Premultiplying a suitable matrix B on both sides of Eq. (34), the reduced row echelon form of matrix P (matrix BP ) can be computed Most C i are still linear combinations of C r i , but C r 14 , C r 15 and C r 17 can be obtained directly. This is because that some rows in the matrix BP have only one none-zero element. The related C r i does not linear combine with the others.
(iv) All remaining about 6.8 × 10 4 sets of L r i give 6.8 × 10 4 sets of C i , but a lot of them give bad convergences. A typical three-flavour ChPT correction at NLO, NNLO and N 3 LO is ∼ 25%, ∼ 7% and ∼ 1.5% [9], respectively. For an observable X, except a 1/2 0 , a 3/2 0 and l r i (i = 2, 3), the following constraints is introduced where the denominator X is their theoretical estimates in Eq. (21). Because the typical correction in each order is only a rough estimate, the upper bounds in Eq. (36) are a slightly more than them.
where both denominators a are their theoretical estimates. The constraints at NLO and the higher order are the same as Eq. (36). For a 3/2 0 , as discussion in Sec. V A, the sign of the LO contribution is opposite to the NLO and the NNLO ones. Its absolute value of the theoretical estimate is much smaller than the one at the LO. Hence, the constraint for a 3/2 0 at NNLO is looser than the other one.
(v) All the remaining C i distribute widely. However, in ChPT, the absolute values of C i can not be very large. A lot of references have estimated the values of C r i (see Table VIII in Appendix B). These values constrain the ranges of C i . Some results which are quite different from those in the literature are excluded. The ranges of C i are chosen as whereC r i are the mean value of C r i in Table VIII and σ C r i are their standard deviations. Some outliers are removed in the calculation. These values can be found in Table IX. We choose the intervals are 5σ C r i wide, because 3σ C r i wide intervals can not give a large enough sets of C r i in our method (see the discussion below). The constraints for m 2 α (α = π, K, η) and F α (α = π, K) in Eqs. (25) and (26) are replaced by the following constraints to constrain the ranges of C r i , where F α /F 0 is the theoretical estimate of the decay constants. Not all C r i are satisfied these constraints. The corresponding L r i are also removed.
(vi) The distributions of most remaining C i are similar to normal distribution, only a few of them have a little asymmetry. The mean values and standard deviations of C i are regarded as their estimates and errors, respectively. To this step, 13114 sets of L r i and C i are left. This number is large enough in statistics. To save computation time, we only use the mean values of C i to estimate C r i , except C r 14 , C r 15 and C r 17 (they have been obtained in Eq. (35)). Eqs. (34) and (38) constrain a high-dimension area in the space of C r i . The centre of mass for this area is regarded as the estimated value of C r i , the standard deviation of this area is regarded as the error of C r i . We use Monte-Carlo method to determine them. First, the subspace satisfied Eq. (34) can be parameterized by 35 − 17 = 18 parameters. The upper and the lower boundaries in Eq. (38) limit this subspace to an 18 dimensional convex polyhedron. Then one can generate numerous random points in a 18 dimensional parallel polyhedron by computer. The parallel polyhedron needs large enough in order to cover the convex polyhedron. Finally, the points in the convex polyhedron are picked out. The mean value of these points is the estimate of C r i , and the standard deviation of these points is the estimate of ∆C r i . We find that 3σ C r i wide intervals are too narrow to generate points, so 5σ C r i wide intervals are chosen.
(vii) Using the values of C i , L r i can be determined by Method I.

VII. THE RESULTS BY METHOD II
The estimates and errors of C i are given in the Table IV. If we randomly select a half of C i , the mean values and standard deviations are unchanged. It indicates that the number of the samples is sufficient. Most of their relative errors are small enough, only ∆ C i / C i (i = 1, 2, 3, 10) have large values. The distributions of C r i are shown in Fig. 1 and Fig. 2 To obtain a more precise value, one needs more reasonable constraints or more inputs. Eleven C r i (i = 11, . . . , 16,19,21,25,31,34) are nearly boundary independent. They are more reasonable. C r i (i = 11, 12, 13) are very especial. There exists a cliff in the centre. We solve the linear programming problem (Eqs. (34) and (38)) and find that there exists no solution in the other side of the cliff. This is because the other boundaries of C r i limit the solution. However, the cliff is far away from both boundaries. These LECs are also boundary independent.
The results of C r i are shown in Table V. Actually, in the calculation, the linear combinations C r 63 − C r 83 + C r 88 /2 and C r 66 − C r 69 − C r 88 + C r 90 arise as a whole. We also present them in Table V. Our initial ranges of C r i are more wide and many sets of C r i can give a small value of χ 2 /dof. Hence, their errors seem large. Some results in Ref. [9] mark with an asterisk, these values are very closed to the original data in the website [21] (less than 10 −10 ). We guess these results meet the fitting boundaries. The symbol "≡ 0" in Column 4 and Column 8 means these values are zeros in the large-N c limit. However, most intervals of C r i (i = 2, 6, 11, 13, 15, 23) are far away from zero. It indicates that these LECs do not satisfy the large-N c limit, that is why the results in Ref. [16] can not fit L r i well by Method I. The second column in Table VI lists the results for L r i . Comparing with Ref. [9], χ 2 /dof = 4.3/9 is closer to 1. It means that the constraints in Sec. VI relieve the overfitting problem. Convergences play an important role in the fit. For the normal distribution in Eq. (32), the initial standard deviation of L r i = 11.5%. Now, most relative deviations (in the third column) are less than 11.5% obviously, except for L r 2,3 . Their values are a little large. We consider they   The results of C r i in units of 10 −6 . The brackets "[" and "]" mean the results are dependent on the lower and the upper boundaries, respectively. The parentheses "(" and ")" mean the results are independent on the lower and the upper boundaries, respectively. The results with an asterisk mean the original data in Ref. [21] is very closed to those in Ref. [9] (less than 10 −10 ). The symbol "≡ 0" for the results in Ref. [16] means these values are zeros in the large-Nc limits.

LECs results
Ref. [9] Ref. [16] LECs results Ref. [9] Ref. [ are also acceptable. Pct L r 2 and Pct L r 3 in the third column in Table VI are larger than the others. The main reason is that a set of L r i containing small L r 2 (large L r 3 ) are much easier to be picked out in the Step ii (iv). For comparison, Column 4 presents the average values of the obtained sets in Step vi in Sec. VI, Column 5 lists the relative deviations between Column 4 and Column 6. The values in Column 2 and Column 4 are closed to each other. It seems that averaging C r i or averaging L r i first is nearly no difference. These results are also not much different from Fit 2 and the results in Ref. [9]. Because the results in Column 2 are related to C r i in Table V, we choose them as our L r i results.  In Table VII, the values of the observables at each order and Pct order are listed. It can be seen that most observables have a good convergence, except a 1/2 0 and a 3/2 0 at NNLO. They have a poor convergence, because the contributions at NNLO need large values (see the Step iv in Sec. VI). Whether the higher-order values are really small or not requires a more reasonable analysis. It is beyond this work. The truncation errors in the fifth column are very small, all Pct HO are less than 4%, except l r 2 . However, the absolute value of l r 2 decreases order by order. It is not a contradiction. Now, with the hypothesises in Sec. II, the first three problems in the introduction are solved, and the cause of the last two problems is also found.

VIII. SUMMARY
In this paper, we have computed the NLO and the NNLO LECs for pseudoscalar mesons with a new method (Method II). The results are present in Table VI and Table V, respectively. The higher-order contributions are considered in the computation with some hypotheses, i.e. the theoretical values of observables are satisfied the convergence in ChPT, all L r i are stable and C r i are consistent with those in the other references. The results nearly satisfy all these hypotheses. All random processes are repeated several times. The results are nearly unchanged. They are reasonable. Some linear combinations of C r i called C r i are also given. Their values are more reliable. If one knows the more exact values of some C r i , some other C r i can be obtained by these C r i . First, a modified global fit method is used to obtain L r i . If they are only fitted at NLO, the results are very closed to the NNLO fitting results in Ref. [9]. It indicates that the higher-order estimates have a good prediction. They are reliable. However, some L r i deviate from the NLO fitting ones in Ref. [9] too much. The main reason is that the higher-order estimates of f s , g p , a 0 0 and a can not be fitted well, because their LO contributions can not give a good prediction and the NLO contributions tend to be small. For the NLO fit, all the theoretical values of observables have a good convergence. However, at NNLO, we have tried two sets of C r i in the references, but the results are not very good. Later, we use a new method to obtained both L r i and C r i . The idea is that the linear independent combinations of C r i ( C i ) are obtained first, and then C r i are estimated with Monte-Carlo method, finally L r i are fitted by these C r i . Some C r i are dependent on the initial boundaries. In order to obtain more precise values, these C r i need more information to narrow the boundaries. The other C r i are boundary-independent. They can be limited in some reliable intervals. The relative errors between the NLO fitting results and the one by this method are all small. All observables have a good convergence, except a The constraints in this paper are very weak. Hence, some error bars of the NNLO LECs are large. If another method can introduce more restrict constraints, their error bars may be narrower. We hope that this new method not only determines the LECs in ChPT for mesons, but it will also generalize to ChPT for baryons and another effective field theory in the future.