Constraints on low energy QCD parameters from $\eta \to 3\pi$ and $\pi\pi$ scattering

The $\eta \to 3\pi$ decays are a valuable source of information on low energy QCD. Yet they were not used for an extraction of the three flavor chiral symmetry breaking order parameters until now. We use a Bayesian approach in the framework of resummed chiral perturbation theory to obtain constraints on the quark condensate and pseudoscalar decay constant in the chiral limit. We compare our results with recent CHPT and lattice QCD fits and find some tension, as the $\eta \to 3\pi$ data seem to prefer a larger ratio of the chiral order parameters. The results also disfavor a very large value of the pseudoscalar decay constant in the chiral limit, which was found by some recent works. In addition, we present results of a combined analysis including $\eta \to 3\pi$ decays and $\pi\pi$ scattering and though the picture does not changed appreciably, we find some tension between the data we use. We also try to extract information on the mass difference of the light quarks, but the uncertainties prove to be large.


Introduction
Spontaneous breaking of chiral symmetry (SBχS) is a prominent feature of the QCD vacuum and thus its character has been under discussion for a long time [1,2]. The principal order parameters are the quark condensate and the pseudoscalar decay constant in the chiral limit 1 2) where N f is the number of quark flavors q considered light and m q collectively denotes their masses. A a µ are the QCD axial vector currents, while F a P the decay constants of the light pseudoscalar mesons P. The two flavor parameters are usually denoted as Σ and F, while the three flavor ones as Σ 0 and F 0 .
Chiral perturbation theory (χPT) [3,4,5] is constructed as a general low energy parametrization of QCD based on its symmetries and the discussed order parameters appear at the lowest order of the chiral expansion as low energy constants (LECs). Interactions of the light pseudoscalar meson octet, the pseudo-Goldstone bosons of the broken symmetry, directly depend on the pattern of SBχS and thus can provide information about the values of Σ(N f ) and F(N f ).
A convenient reparametrization of these order parameters, relating them to physical quantities connected with pion two point Green functions, can be introduced [2] Z( Z(N f ) = 0 would correspond to a restoration of chiral symmetry and X (N f ) = 0 to a case with vanishing chiral condensate. Standard approach to chiral perturbation series tacitly assumes values of X (N f ) and Z(N f ) not much smaller than one, which means that the leading order terms should dominate the expansion. Several recent results for the two and three flavor order parameters are listed in Tables 1 and 2, respectively. As can be seen, while the two flavor case is quite settled, the values of X (2) and Z(2) indeed being not much smaller than one, the situation in the three flavor sector is much less clear. Some analyses suggest a significant suppression of X(3) and/or Z(3) and thus a non-standard behavior of the spontaneously broken QCD vacuum. It can also be noted that the lattice averaging group FLAG [6] does not report an average for the three flavor chiral order parameters.
Up, down and strange quark masses are other parameters with strong influence on low energy QCD physics. A commonly used reparametrization can be introduced (1.4) 1 We will abbreviate these to chiral condensate and chiral decay constant in the following The values for the light quark mass average and the strange to light quark mass ratio are well known from lattice QCD and QCD sum rules [6,15]. On the other hand, as can be seen in Table 3, the isospin breaking parameter R, directly related to the light quark mass difference, has not been determined with such a high precision by these or other methods yet.
In this paper, we use a Bayesian approach in the framework of resummed chiral perturbation theory to extract information on the three flavor chiral condensate, chiral decay constant and the mass difference of the light quarks. Our experimental input are well known observables connected to η → 3π decays and ππ scattering. We assume a reasonable convergence of Green functions connected to these observables and investigate the constraints this assumption can provide for the discussed parameters. The results presented here are a significant update on our initial reports [20,21,22] In Section 2 we shortly summarize our theoretical foundation. Section 3 discusses the η → 3π decays, while Section 4 provides an overview of our calculation of these processes. Section 5 introduces ππ scattering into our analysis. The Bayesian statistical approach is reviewed in Section 6 and a discussion of our assumptions can be found in Section 7. Section 8 is concerned with investigating the compatibility of our theoretical predictions with the ππ scattering data. We employ a χ 2 based analysis in Section 9 to evaluate the quality with which our theoretical predictions reconstruct the experimental data and use it to choose between several assumptions. Finally, in Section 10, the main results of the Bayesian analysis are presented and compared with available literature. We conclude in Section 11. phenomenology Z(2) X(2) ππ scattering [7] 0 NNLO χPT (BE14) [9] 0.59 0.63 NNLO χPT (free fit) [9] 0.48 0.45 NNLO χPT ("fit 10") [10] 0.89 0.66 ReχPT ππ+πK [11] >0   [16] 37 η → 3π NNLO χPT [16] 41.3 η → 3π dispersive [17] 37.7±2.2 η → 3π dispersive [18] 34.2±2.2 η → 3π dispersive [19] 32.7±3.0 35.7±2.6 Table 3: Chosen results for the isospin breaking parameter R. r = 37.3 ± 0.34 [6] was used to obtain R from Q in ref. [18] and [19].

Resummed χPT
We use an alternative approach to chiral perturbation theory, dubbed resummed χPT (ReχPT) [23], which was developed in order to accommodate the possibility of irregular convergence of the chiral expansion. This is a typical scenario if the X (3) and Z(3) were indeed suppressed and the leading order was not dominant in the chiral expansion. In such a case one would have to be careful in the way how chiral expansion is defined and dealt with, as reshuffling of chiral orders could lead to unexpectedly large higher orders.
The procedure can be shortly summarized in the following way: • The Standard χPT Lagrangian and power counting [3,4,5] is used. In particular, the quark masses m q are counted as O(p 2 ). • Only expansions of quantities related linearly to Green functions of QCD currents are trusted, these are called safe observables. It is assumed that the next-to-next-to-leading and higher orders in these expansions are reasonably small, though not necessary negligible. Leading order terms do not need to be dominant. • Calculations are performed explicitly to next-to-leading order, higher orders are included implicitly in remainders. The first step consists of performing the strict chiral expansion of the safe observables, which is understood as an expansion constructed in terms of the parameters of the chiral Lagrangian, while strictly respecting the chiral orders. • In the next step, the strict expansion is modified in order to correct the location of the branching points of the non-analytical part of the amplitudes, which need to be placed in their physical position. This can be done either by means of a matching with a dispersive representation or by hand. The bare expansion is obtained. • After that, an algebraically exact non-perturbative reparametrization of the bare expansion is performed. It is obtained by expressing the O(p 4 ) LECs L i in terms of physical values of experimentally well established safe observables -the pseudoscalar decay constants and masses. The procedure generates additional higher order remainders. We refer to these as indirect remainders.
• The physical amplitude and other relevant observables are then obtained as algebraically exact non-perturbative expressions in terms of the related safe observables and higher order remainders. • The higher order remainders are not neglected, but estimated and treated as sources of error.
The hope for resummed χPT is that by carefully avoiding dangerous manipulations a better converging series can be obtained. The procedure also avoids the hard to control NNLO LECs by trading them for remainders with known chiral order.
In general, a chiral expansion of a safe observable is written in the following way: G is the remainder which contains all terms starting with NNLO. The basic assumption of this paper is that δ (6) G ≪ 1 for our chosen observables. We will quantify this assumption later on.

η → 3π decays
The η → 3π isospin breaking decays have not been exploited for an extraction of the chiral order parameters so far, yet we argue there is valuable information to be had. The theory seems to converge slowly for the decays. One loop corrections were found to be very sizable [24], the result for the decay width of the charged channel was 160±50 eV, compared to the current algebra prediction of 66 eV. However, the experimental value is still much larger, the current PDG value is [25] Γ + exp = 300 ± 12 eV. (3.1) The latest experimental value for the neutral decay width is [25] Γ 0 exp = 428 ± 17 eV. (3.2) Only the two loop χPT calculation [16] has succeeded to obtain a reasonable result for the widths.
As we have shown in [26], we argue that the four point Green functions relevant for the η → 3π amplitude (see (4.1) below) do not necessarily have large contributions beyond next-to-leading order and a reasonably small higher order remainder is not in contradiction with huge corrections to the decay widths. The widths do not seem to be sensitive to the details of the Dalitz plot distribution, but rather to the value of leading order parameters -the chiral decay constant, the chiral condensate and the difference of u and d quark masses, i.e. the magnitude of explicit isospin breaking. Moreover, access to the values of these quantities is not screened by EM effects, as it was shown that the electromagnetic corrections up to NLO are very small [27,28]. This is our motivation for our effort to extract information about the character of the QCD vacuum from this decay.
The Dalitz plot distributions are experimentally well known as well [29,30,31,32,33]. The usual parametrization of the square of the amplitude is defined as |A(s,t; u)| 2 = |A(s 0 , s 0 ; s 0 )| 2 1 + ay + by 2 + dx 2 + f y 3 + gx 2 y + . . . (3.3) in the charged decay channel and as |A(s,t; u)| 2 = |A(s 0 , s 0 ; s 0 )| 2 (1 + 2αz + . . .) (3.4) in the η → 3π 0 case, where and s 0 = 1 3 (s + t + u) is the center of the Dalitz plot. However, as we have discussed in detail in [26], we have not found the convergence of the theory in the case of the slopes reliable enough to include all the experimentally measured Dalitz plot parameters into the analysis. To stay on the conservative side, we used the lowest order parameter a in the charged channel only. The latest and most precise experimental value is [33] a = −1.095 ± 0.004. (3.6)

Calculation
The details of the calculations with explicit formulas can be found in [26]. We only summarize the basic steps here, which closely follow the procedure outlined in [34].
We start by expressing the charged decay amplitude in terms of 4-point Green functions G i jkl , obtained from the generating functional of the QCD currents. We compute at first order in isospin breaking, the amplitude then takes the form where ∆ G D is the direct higher order remainder to the 4-point Green functions. The physical mixing angles to all chiral orders and first order in isospin breaking can be expressed in terms of quadratic mixing terms of the generating functional to NLO (M In this approximation the neutral decay amplitude can be straightforwardly obtained from the charged one using isospin symmetry and charge conjugation invariance In accord with the method, O(p 2 ) parameters appear inside loops in the strict chiral expansion, while physical quantities in outer legs. Such a strictly derived amplitude has an incorrect analytical structure due to the leading order masses in loops, cuts and poles being in unphysical positions. We have developed several ways to account for this in [34], explicit formulas for the η → 3π case are listed in [26]. The simplest approach is to exchange the LO masses in unitarity corrections and chiral logarithms for physical ones. A more sophisticated method is to use a dispersive representation for the unitarity part of the amplitude where U (s,t; u) is the unitary part and P(s,t; u) is an unknown polynomial. This can be obtained by using the reconstruction theorem [35], first used in [36]. Then the two representations can be sewed together where the polynomial part G pol (s,t; u) is obtained from (4.1) by the sewing procedure. However, there is an ambiguity in the derivation of U (s,t; u). When using the Cutkosky rule, there is a freedom in the way how to define the relation of the amplitude and the Green function of the entering sub-process at leading order. The most straightforward way is to keep the order by order relation Such a definition has the advantage of satisfying perturbative unitarity. On the other hand, a suppression factor F 4 0 /(∏ i F 2 P i ) appears in the loop functions, where P i are the pseudocalars running in the loop. Because we expect the unitarity correction in the case of the η → 3π decays to be sizable, we decided to implement an alternative definition No suppression factor occurs in this case. Both approaches are valid and they differ in the definition of the higher order remainder ∆ G D . In our view, we should prefer such a definition where the higher orders are under better control. The approaches will be numerically compared in Section 9. The next step is the treatment of the LECs. As discussed above, the leading order ones, as well as quark masses, are expressed in terms of convenient parameters At next-to-leading order, the LECs L 4 -L 8 are algebraically reparametrized in terms of pseudoscalar masses, decay constants and the free parameters X , Z and r using chiral expansions of two point Green functions, similarly to [23]. Because expansions are formally not truncated, each generates an unknown higher order remainder We don't have a similar procedure ready for L 1 -L 3 at this point. Therefore we collect a set of standard χPT fits [9,10,37,38] and by taking their mean and spread, while ignoring the much smaller reported error bars, we obtain an estimate of their influence These uncertainties enter our statistical analysis. However, as it is discussed in [26], the results depend on the value of the constants L 1 -L 3 only very weakly.
The O(p 6 ) and higher order LECs, notorious for their abundance, are implicit in a relatively smaller number of higher order remainders. We have eight indirect remainders -three generated by the expansions of the pseudoscalar masses, three by the decay constants and two by the mixing angles. We expand the direct remainder to the 4-point Green functions around the center of the Dalitz plot s 0 = 1/3(M 2 η +2M 2 π + +M 2 π 0 ) to second order in Mandelstam variables and thus get four derived direct remainders, two NLO ones (∆ C , ∆ D ) and two NNLO ones (∆ A , ∆ B ). We should note that in this approximation we completely miss the analytical structure of the amplitude at higher chiral orders, which includes ππ rescattering effects at NNLO and higher. A deeper discussion of the issue can be found in [26]. We argue that because the experimental curvature of the Dalitz plot is very small [29], the expansion to second order in the Mandelstam variables is sufficient for the purpose of calculating the decay widths and the lowest order Dalitz slope a. On the other hand, the theory seems unreliable for the higher order Dalitz parameters, especially in the case of b and α, thus we have chosen to avoid them in this analysis. For the calculation of the decay widths, we need to numerically integrate over the kinematic phase space. In order to perform the computation efficiently enough, we are forced to expand the amplitude around the center of the Dalitz plot (4.14) The same argument as above hold here as well -the curvature of the Dalitz plot is tiny, therefore this is a very good approximation for the objective of calculating the decay widths.

ππ scattering
In addition to the η → 3π parameters discussed above, we employ ππ scattering in a very similar way to [23]. We use the two lowest order subthreshold parameters in the expansion of the polynomial part of the amplitude, α ππ and β ππ is the unitary part of the amplitude to NNLO, not given here explicitly. α ππ and β ππ can be expressed as We use the experimental values extracted in [7] α exp ππ = 1.381 ± 0.242, β exp ππ = 1.081 ± 0.023, ρ ππ = −0.14. ρ ππ is the correlation coefficient between the two parameters.

Bayesian statistical analysis
For the statistical analysis, we use an approach based on Bayes' theorem [23] P( where P(X i |data) is the probability density of the parameters and remainders, denoted as X i , having a specific value given the observed experimental data.
In the case of experimentally independent observables, P(data|X i ) is the known probability density of obtaining the observed values of the included observables O k in a set of experiments with uncertainties σ k under the assumption that the true values of X i are known Our observables are the charged and neutral decay widths and the Dalitz slope a of the η → 3π decays and the ππ scattering subthreshold parameters α ππ and β ππ . However, the pair of ππ scattering observables cannot be treated as independent and we need to introduce a correlated probability function for this sector where (6.4) We resort to Monte Carlo sampling in order to perform the numerical integration in (6.1). It turned out that the uncertainty of latest experimental measurement of the Dalitz parameter a (3.6) is so small that it is in fact a complication for performing the numerical integration. From this point of view the experimental error is negligible and therefore we can model the experimental distribution as a δ function rather then a normal distribution The experimental data thus become a constraint which can be solved algebraically. The solution is most straightforward in the case of the remainder δ B where The remainder δ B becomes fixed by this constraint and is thus no longer a source of uncertainty. P(X i ) in (6.1) are the prior probability distributions of X i . We use them to implement the theoretical uncertainties connected with our parameters and remainders. In this way we keep the theoretical assumptions explicit and under control. It also allows us to test various assumptions and formulate if-then statements as well as implement additional constraints (see below).

Assumptions
The following list summarizes the higher order remainders we need to deal with: We use an estimate based on general arguments about the convergence of the chiral series [23] δ (4) G collectively denotes all the remainders listed above, with the exception of δ B , which is fixed by the constraint (6.6) in the main analysis. We implement (7.1) by using a normal distribution with µ = 0 and σ = 0.3 or σ = 0.1 for the NLO or NNLO remainders, respectively. The remainders are thus limited only statistically, not by any upper bound.
The remainders δ α ππ and δ β ππ , connected with ππ scattering, are defined in such a way that they are of order O(mm s ) instead of O(m 2 s ). This follows from the particular form of the chiral expansion of the subthreshold parameters α ππ and β ππ in them → 0 limit, see [23] for details. Therefore we can consider these remainders to be of the order δ α ππ ≈ ±0.03, δ β ππ ≈ ±0.03.
We assume the strange to light quark ratio r to be known and use the lattice QCD average [39] At last, we are left with three free parameters: These control the scenario of spontaneous breaking of chiral symmetry and isospin breaking in our results. In the case of X and Z, we use a constraint from the so-called paramagnetic inequality [2] and assume these parameters to be in the range For the two flavor order parameters, we use the lattice QCD values [8]. In addition, similarly to [23], we implement constraints following from We use two approaches to deal with R. In the first one we assume it to be a known quantity. We use the N f =2+1 lattice QCD average [39] R = 35.8 ± 2.6. (7.7) Alternatively, we leave R free, or more precisely, assume it to be in a wide range R ∈ (0, 80).

Subthreshold parameters of ππ scattering
As mentioned in Section 3, in [26] we tested the compatibility of a reasonable convergence of the chiral series, laid out explicitly in the form of assumptions in the previous section, with the experimental data in the case of the η → 3π observables. This is an important first step in order to avoid using observables which are problematic from the theoretical or experimental point of view. Analogously, in this section we take a closer look at the subthreshold parameters α ππ and β ππ , which was not done in [23].
We numerically generated a large number of theoretical prediction for α ππ and β ππ , statistically distributed according to the assumptions described in Section 7. The parameter Z was fixed in two scenarios (Z=0.5 and Z=0.9), while X was varied in the full range 0 < X < 1. Figure 1 displays the obtained theoretical distributions in comparison with the experimental data.
As can be seen, both parameters show a significant dependence on X , while β ππ depends on Z as well. In the case of β ππ , a broad range of the generated theoretical predictions is consistent with experimental data. Even though the observable seems sensitive to the values of the chiral order parameters, the amount of available information to be extracted might be limited by the experimental error, which is quite substantial.
The picture seems to be more tricky when considering α ππ . The central experimental value is outside the 2σ band of the theoretical distribution in the whole range of the free parameters. If taken at face value, possibly confirmed by more precise data, it would indicate a very small value of X , i.e. a vanishing chiral condensate. Such a scenario would be, however, inconsistent with any current determination of the order parameter (see Tab.2). 2 The experimental error on the value of α is very large and our prior expectation therefore is that a large correction of the central value is very possible. It would be thus advisable to be cautious when interpreting the outcomes based on this data in the following.
Our conclusion with regard to the suitability of the subthreshold parameters of ππ scattering for the purpose of extraction of information about the chiral order parameters is hence twofold -both α ππ and β ππ seem sensitive to the value of one or both of them, but at present available information is limited due to the quality od experimental data we have at hand.

χ 2 based analysis
In addition to the Bayesian analysis described in Section 6, the results of which will follow in Section 10, we also perform a search for the minimum of the χ 2 distribution in the Monte Carlo generated set of data points. The aim is to check the quality with which this set of theoretical predictions can reconstruct the experimental data, which is hard to quantify using the Bayesian method. This test also enables us to compare various theoretical approaches, e.g. the alternative ways the dispersive representation can be implemented.
Because we have much more free parameters than experimental inputs, we generally expect a well working theory to be able to reconstruct the experimental data very precisely with a large enough sample of generated data points. This means the minimum of χ 2 /n, where n is the number of experimental observables employed, should be close to zero. On the other hand, a min.χ 2 /n ≈ 1 means there is typically a 1σ deviation between the best of the generated theoretical predictions and the experimental data, which might signal some tension and would not be fully satisfactory.
As the minimum of the χ 2 distribution is subject to fluctuations stemming from the statistical nature of our procedure, we also report the number of points for which χ 2 /n < 1. This metric then reveals how well the region of the parameter space where the experimental data lie is covered by the generated theoretical predictions. A reasonably high number should be considered a necessary condition for a well founded analysis, as it demonstrates that the theory has no obvious problem to reconstruct the experimental data and also that we have generated a large enough sample of points.
When comparing different theoretical approaches, one model can overlay the experimental region with fewer points than the other for several reason. It might be that the points originate from a less probable part of the theoretical distribution (e.g., the tail) which means that less likely values of the parameters and remainders are needed in order to reconstruct the experimental data. In that sense, the model is less probable to be true. The χ 2 based analysis can thus form a basis for preference when choosing between alternative approaches. However, it is also possible that one model is simply more sensitive to the values of some of the parameters or remainders, which means the theoretical distribution has a larger spread. This is no reason to exclude to model, of course, and one should be aware of this possibility and check for it.
Let us first explore the issue of deciding between the two dispersive representations used to define the bare expansion of the η → 3π amplitude, as discussed in Section 4. For this purpose we numerically generated a set of 4 − 6 · 10 6 theoretical predictions and constructed a χ 2 distribution. In order to check the consistency of the generated predictions and a complete set of η → 3π data, we have not used the constraint (6.6), but rather an experimental value a = −1.09 ± 0.02. (9.1) This is a slightly older measurement [29] with large enough uncertainty for our purpose.
free parameters exp.data disp.aproach min.χ 2 /n N(χ 2 /n < 1)  As can be seen in Table 4, the approach (4.7) is able to reconstruct the data more precisely, while the theoretical distributions have a very similar form. This conforms to our intuition from Section 4 and thus, in what follows, we use the representation based on (4.7) exclusively.
We use the PDG [41] as the source of input for the values of pseudoscalar masses and decay constants. Because the isospin breaking parameter 1/R carries the isospin symmetry breaking, we need to choose a value of these constants in the isospin limit. We use the averaged kaon mass and the well known decay constants of the charged pion and kaon. However, we found the situation to be quite subtle for the case of the pion mass. For the charged decay η → π + π − π 0 it seems appropriate to use the averaged mass while in the case of the neutral channel η → 3π 0 to use the neutral pion mass M π 0 . However, in this approach the isposin relation (4.3) is not exactly fulfilled. Alternatively, one could use the same pion mass for both channels, either the averaged or neutral one, and satisfy the relation (4.3). One could argue that the difference in the results for the decay widths should be very subtle, of the order 1/R, and that certainly seems to be true. However, the ratio of the decay widths is known very precisely [25] r Γ = 1.43 ± 0.02. (9.3) This is an indication that even a slight difference in prediction of order 1/R might actually influence whether one can obtain an accurate enough prediction for both the decay widths at the same time.
As we have found out, this is really the case in the approach with an identical pion mass in both the charged and neutral channel amplitudes, where the predicted ratio r Γ comes out too high. This is reflected in both the minimum of χ 2 /n not being quite close to zero and in the number of points for which χ 2 /n < 1 to be substantially lower in our χ 2 based test for any value of the parameters in the allowed range, as can be seen in Table 5. Meanwhile, the form of the distributions do not change considerably, as might be expected when only the numerical value of an input parameters is changed. free parameters exp.data pion mass min.χ 2 /n N(χ 2 /n < 1)  In other words, in this case the theory seems to have a harder time to reproduce the experimental data with the required precision. This might look surprising given the number of free parameters in the fit, but one has to realize that while the decay widths in the two channels are independent from the experimental point of view, they are very strongly correlated on the theoretical side given the relation (4.3). The experiment in fact shows that this relation is not precisely fulfilled in nature. This is the reason we have chosen to implement distinct values of the pion mass in the two channels for the main analysis, which is the same approach as was taken in [16]. Table 6 shows a summary of the χ 2 based tests for our main analysis, as presented in the next section. This uses the dispersive representation (4.7), distinct pion mass values for the two η → 3π decay modes and the constraint (6.6). The total number of generated data points is 2 · 10 7 . free parameters exp.data min.χ 2 /n N(χ 2 /n < 1) Table 6: χ 2 based comparison for the main analysis used in Section 10.
As can be seen, the number of points seems sufficient. The coverage is better in the approach with η → 3π observables only. We can observe a drop in precision when α ππ is included, which correlates with the discussion in Section 8. As we will see, this corresponds to the fact that the signal in the used experimental value of α ππ is in some tension with the one contained in the η → 3π data.

Results and discussion
In this section we present the outputs of the Bayesian analysis, i.e. the probability density functions P(X i |data), where X i ∈ {X , Z,Y, R} are the chiral symmetry breaking and explicit isospin breaking parameters, respectively and data represent a subset of the set {Γ + , Γ 0 , a, α ππ , β ππ } of the η → 3π and ππ → ππ observables.
The calculated probability density distributions P(X i |data) are depicted in Fig.2, Fig.6 and Fig.7. The colors in these figures highlight the confidence regions of the parameter space with a particular confidence level 3 . The corresponding one dimensional probability density functions for the parameters X , Z, Y = X /Z and R are obtained by integrating out all other free parameters. In particular, The results are shown in Fig.3, Fig.4, Fig.5 and Fig.8, along with the priors following from the assumptions (7.5) and (7.6). Error bars in these figures indicate an estimated error of the Monte [13] [14] [9]  Carlo integration, which is sufficiently low in all cases. We summarize some characteristics of these distributions in Tab.7 and Tab.8. Let us first discuss the case with the fixed value of R = 35.8 ± 2.6, which is the lattice QCD average [39] (see Section 7) and include only the η → 3π observables into the analysis. This is what we consider our main set of results.
As can be seen in the left panel of Fig.2, a large part of the parameter space can be excluded at 2σ CL and the parameters X and Z appear to be quite strongly correlated. The latter statement can be made more quantitative by extracting the result for the ratio of the chiral order parameters Y = X /Z = 2mB 0 /M 2 π , see Fig.3 and Tab.7, for which we get: We also obtain a higher bound for the three flavor chiral decay constant   As can be confirmed from Fig.3, a combination of two factors contribute to this result -the η → 3π data disfavor high values of Z and the prior (7.5), induced by the paramagnetic inequality, causes a sharp cutoff at Z = Z(2) = 0.86 ± 0.01. As can be seen in Figure 2, when assuming R = 35.8 ± 2.6, there is some tension with several of the previous determination of the chiral order parameters ( Table 2). The η → 3π data seem to prefer a larger value for the ratio of the order parameters Y = X /Z than recent χPT and lattice QCD fits. In addition, very large values of the chiral decay constant are excluded at 2σ CL. and a relatively small value is favored. The uncertainties, however, are quite large.
RBC/UKQCD+large N c [13] 0.91±0.08 -0.05±0.22 NNLO χPT (BE14) [9] 0.59 0.3 NNLO χPT (free fit) [9] 0.48 0.76±0.18 RBC/UKQCD+ReχPT [12] 0.54±0.06 1.06±0.29 Table 9: Correlation of Z and L 4 in available determinations of F 0 .  Of course, these determinations are hardly compatible among themselves either and provide a very unclear picture. In our view, a reasonable guess could be that there are important differences in how NLO and NNLO low energy constants are treated in various works. In particular, one possible culprit could be the large N c suppressed LEC L 4 , which is known to be anti-correlated with F 0 for a long time [42]. Indeed, this anti-correlation can be observed in the results we quote, see Tab.9.
In our approach, L 4 is reparametrized in terms of the remainders and free parameters, including Z. It can thus vary in a wide range, as we have shown in [26]. As the η → 3π data seem to constrain Z only mildly, as discussed above, we do not get significant information on L 4 either.
When concerning [9], both the main fit (BE14) and the free fit, which are based on the standard χPT at O p 6 and a large set of the experimental observables, are relatively close to the 2σ contour of our distribution. Note however, that the errors of these fits are not at our disposal. The two fits differ precisely in their treatment of L 4 , as indicated in Tab.9. If we roughly estimated the theoretical uncertainty of these fits as the distance of the corresponding points in the (X , Z) plane, then these fits might look quite compatible with our distribution.
The apparent inconsistency with the result of [12] is intriguing. It uses resummed χPT as well, paired with lattice data. One distinction is a different approach to the remainders -the authors use uniform distribution of the remainders inside a closed interval and thus a sharp cutoff. One could speculate that a normal distribution with unbounded tails, as we use, might lead to larger error bars.
The work [13], which is based on a large N c motivated approximation of the standard O p 6 χPT calculations, used on lattice data, reports a very large value of the chiral decay constant and a very low value of L 4 . This is consistent with the large N c picture assumed in this paper, but quite far away from other determinations and in tension with our limit (10.2). In fact, a large part of the region covered by the fit [13] is excluded by our prior for Z, namely by the constraint stemming from the paramagnetic inequality Z < Z(2) (7.5).
Let us now add the ππ → ππ data into the analysis. As can be seen in the right panel of Fig.2 and Fig.4, the picture does not change appreciably when including the subthreshold parameters of ππ scattering. Though a bit disappointing, this outcome is not unexpected considering the significant errors connected with the experimental values of these observables and the weak constraints obtained in [23] and [11].
There is one difference, however, we can observe a slight shift of our probability distribution towards a lower ratio of order parameters Y = X /Z, as is confirmed by the mean Y = 1.2 (P(Y |Γ + , Γ 0 , a, ππ) in Tab.7). This could be interpreted as a move of our predictions in the direction of better compatibility with some of the available determinations (Tab.2). Here we have to be rather cautious though, because, as we have discussed in Sect.8 and Sect.9, we can expect some tension between the two sets of data. And indeed, this can be demonstrated when one compares the obtained probability distributions and confidence intervals for Y from η → 3π (Fig.3, P(Y |Γ + , Γ 0 , a) in Tab.7) with one from ππ scattering alone (Fig.5, P(Y |α ππ , β ππ ) in Tab.7), which are barely compatible. We can conclude that we need a more precise determination of the value of the ππ → ππ subthreshold parameters, especially for α ππ , to be able to provide a more definite outcome and hopefully solve this puzzle of experimental data pointing in opposite directions.
The results with R left as a free parameter are shown in Fig.6 and Fig.7. The uncertainties are large and thus it's hard to constrain R without additional information on the chiral order parameters and the remainders. Even in this case a part of the parameter space can be excluded at 2σ C.L. though. The obtained value for R (Fig.8, Tab.8) is compatible with available results (Table 3).
We can also evaluate the obtained probability densities for X , Z and Y with R left unconstrained (Tab.8, Fig.8). Note that dismissing the very clear information on R is not a reasonable assumption, we use it only as a test of robustness. And in this case, the results for X and Z seem  On the other hand, as Fig.7 shows, R and the ratio Y = X /Z seem to be quite strongly correlated. This also provides some additional insight into the large value of Y obtained from η → 3π data when fixing R = 35.8 ± 2.6. Furthermore, as the ππ scattering data we use drag Y down to lower values, one can observe a correlated shift in probability densities of R to smaller numbers in Fig.7, Fig.8 and Tab.8.

Conclusions
To summarize, we have used statistical methods in the framework of resummed chiral perturbation theory to generate large sets of theoretical predictions for η → 3π and ππ scattering observables, dependent on a variety of parameters and assumptions, and confronted them with experimental data.
We have developed a χ 2 based analysis, which allowed us to form a basis for preference when choosing between alternative assumptions or models. In particular, it showed us that an approach using different pion masses for the charged and neutral decay channel observables in the isospin limit is significantly more consistent with data, despite violating the isospin relation, than using identical pion mass for both η → 3π decay modes.
For the main analysis, we have used Bayesian inference to obtain constraints on the values of three flavor chiral order parameters -the chiral decay constant F 0 and the chiral condensate Σ 0 , which are connected with the spontaneous breaking of chiral symmetry in QCD.
When fixing the light quark difference by input from lattice QCD and using η → 3π observables only (the decay widths in both channels and the Dalitz plot parameter a), we could exclude a large part of the parameters space at 2σ CL and have observed some correlation between the chiral order parameters. We have obtained an upper bound for the chiral decay constant, F 0 < 81MeV at 2σ CL, and have extracted a fairly large value for the ratio of the order parameters Y = 2mB 0 /M 2 π . We have found some tension with several of the previous determination of the chiral order parameters, which, however, are neither very consistent with each other. The picture remains unclear, possibly stemming from differences in assumptions about the low energy constants, the large N c suppressed LEC L 4 being one candidate.
The picture have not changed appreciably when we performed a combined η → 3π and ππ → ππ analysis. However, we have observed some tension between η → 3π and ππ scattering data, which limited our ability to draw more definite conclusions. The possible source is the experimental error of the observables we used, the subthreshold parameters α ππ and β ππ . Specifically, the value of α ππ proved suspicious in our investigation, as it prefers a very low value for the chiral condensate, which is not very consistent with current expectations.
We have also tried to extract information on the difference of light quark masses, but the uncertainties proved to be very large. The result is consistent with available data though.