text:
\section{Introduction}
The masses of the light quarks are fundamental free parameters of the standard
model.
Since quarks are confined inside hadrons, there is no direct method for their
measurement.
The only method of determination is a comparison of the theoretical
prediction for some observable that depends on the quark masses with the
corresponding experimental value. For that end we need a framework, in which the
quark masses occur explicitly, and which can make predictions for such
observables with sufficient precision. Because of quark confinement and the
fact that these masses are very small in comparison to the typical hadron
scales, perturbative quantum chromodynamics (QCD) cannot play such a role
and we need to employ a non-perturbative method. Nowadays, the natural
candidates for such approaches are lattice QCD \cite{Wilson:1974sk} and chiral
perturbation theory (ChPT) \cite{Weinberg:1978kz,GL-scatlen,Gasser:1984gg}.

While the mass $m_s$ and the isospin average mass of $m_u$ and $m_d$, which is
denoted as
\begin{equation}\label{hatm}
\hat{m}=\frac{m_u+m_d}{2}\,,
\end{equation}
have become accessible through the recent lattice simulations (among others \cite{RBC-UKQCD,MILC,HPQCD,PACS-CS}) and agree well with the independent determination via QCD sum rules \cite{Maltman:2001nx,Narison:2005ny,Dominguez:2008jz}, the extraction of the individual masses $m_u$ and $m_d$ from the lattice is still polluted by various
simplifications of the electromagnetic effects that have to be made in
isospin-breaking simulations (cf.\ e.g.\ \cite{Colangelo:2010et}). Therefore, if
we want to determine the individual masses $m_u$ and $m_d$, at the moment ChPT seems to be
the more promising approach.

The most suitable processes for studies of isospin breaking within ChPT in the mesonic sector are
$\eta\to3\pi$ decays. Since these decays violate G-parity,\footnote{
Or equivalently, the decay would be forbidden as a result of isospin conservation and charge conjugation invariance (C-invariance). Indeed, the final state has to have the total isospin $I=0$ and is therefore totally antisymmetric with respect to the permutation of the three pions (the only allowed state is then $\pi^+\pi^-\pi^0$). Due to Bose symmetry, the corresponding amplitude is then totally antisymmetric under exchanges of the momenta of these three pions. On the other hand, according to C-invariance, the amplitude has to be symmetric with respect to the exchange of the momenta of the $\pi^+$ and the $\pi^-$, which implies that the amplitude is zero.} they have to proceed via isospin-breaking effects. There are two mechanisms of this breaking, either through the
electromagnetic (EM) interactions,
which are proportional to the electric charge squared,
or through the isospin-breaking mass difference between $u$ and $d$ quarks,
\begin{equation}\label{eta hamiltonian QCD}
\mathcal{H}^{IB}_{\mathrm{QCD}}(x)=\frac{m_d-m_u}{2}\left(\bar{d}(x)d(x)-\bar{u}
(x)u(x)\right).
\end{equation}

Even though the EM interactions have a sizable effect on the difference $m_{\pi^\pm}-m_{\pi^0}$
and on the pion decay constant $F_\pi$, it has turned out that their influence
on the $\eta\to3\pi$ decay amplitudes is very small
\cite{Sutherland:1966zz,Baur:1995gc,Ditsche:2008cq}. Hence,
$\mathcal{H}^{IB}_{\mathrm{QCD}}$ represents the dominant contribution and the
amplitude is proportional to $m_u-m_d$, which is usually presented in one of the
following ratios
\begin{equation}\label{RQ}
R = \frac{m_s - \hat m}{m_d-m_u}, \qquad
\frac{1}{Q^2}=\frac{m_d^2-m_u^2}{m_s^2-\hat{m}^2}
\end{equation}
that are connected via ($r=\frac{m_s}{\hat{m}}$)
\begin{equation}\label{RQ relace}
Q^2=\frac{1}{2}R(r+1).
\end{equation}

Consequently, a measurement of the decay rates of the $\eta\to 3\pi$ processes enables us a direct access to this difference (and by the use of the values $\hat{m}$ and $m_s$
from the lattice also to the individual masses of these two lightest quarks).
Of course, in order for this extraction to be possible, it is necessary to have ChPT predictions for these decay
rates with a sufficient degree of accuracy.

Achieving this is, however, a non-trivial task. The tree-level predictions,
which are equivalent to the PCAC results
(e.g.~\cite{Bell:1996mi,Cronin:1967jq,Osborn:1970nn}), would indicate a very
large difference between $m_u$ and $m_d$. Furthermore, the true energy dependence of
the amplitudes is definitely different from the trivial one that PCAC proposes. The
sizable one-loop corrections \cite{Gasser:1984pr} were still not sufficient to
correct these discrepancies. At last, the inclusion of the two-loop corrections
\cite{BG}, which are also sizable, led to the predictions for both $R$ and the Dalitz
parameters describing the energy dependence of the amplitudes
(cf.~Tables~\ref{tabulka Dalitz etac} and \ref{tabulka Dalitz eta0} below) that
were in reasonable agreement with expectations.

Nevertheless, if we study these results in greater detail, we find some hints
that the \mbox{2-loop} ChPT determination of $R$, which we are interested in, can still be
inaccurate. The feature often put forward in this respect is the discrepancy between
the experimentally measured and the predicted values of the Dalitz parameters (defined
in Sec.~\ref{Dalitz plot parametrization}), mainly for the neutral parameter $\alpha$ of (\ref{dalitz eta0}).
For a better quantification of the difference between experiment and theory, let us introduce
\begin{equation}\label{chi2}
 \chi^2 \equiv \Bigl( \frac{\text{exp} - \text{theory}}{\sigma[\text{exp}]}
\Bigr)^2\,,
\end{equation}
where the theory enters only via its central value.
Using this quantity when comparing the prediction of \cite{BG}, $\alpha = 0.013$, with
the best measurement of this observable by MAMI-C \cite{MAMI-C}, $\alpha = -0.032\pm
0.003$, we obtain indeed a huge difference of $\chi^2 = 225$.
However, there is a parameter for which this discrepancy is even more apparent, namely
$b$ from (\ref{dalitz etac}). Comparing the ChPT value $b=0.394$ with the measurement by KLOE \cite{KLOE-etac} $b=0.124\pm0.012$ produces
$\chi^2\sim 500$. This raises the question about the origin of these discrepancies, and whether
and to which extent they can also affect the determination of $R$.

As was already stressed in \cite{BG}, the explanation of this difference between theory and experiment
 can be provided by the large theoretical error bars presented there (thereby making the theoretical and the experimental values compatible).
The non-renormalizability of ChPT represents a major drawback of this theoretical framework when it becomes necessary to include higher and higher orders. Indeed, including two-loop effects to the computation means a rapid
increase\footnote{Note that $O(p^6)$ ChPT contains 102 (2+10+90) free
LECs. In order to make this theory at the given order predictive, we would thus
need to make at least 102 measurements determining these constants. Obviously,
not all of these constants appear in a given amplitude --- only a
subset of them contributes to $\eta\to3\pi$, see below (Sec.~\ref{sec:chpt}).} of the number of \emph{a priori}
unknown low energy constants (LECs) that have to be estimated before we can
get a reliable prediction. We are far from a determination of all required
LECs from experiment (or lattice), and hence for many of them we have to rely
on some estimates, predominantly of the resonance saturation type
\cite{Ecker:1988te,Ecker:1989yg,Moussallam:1994at,Knecht:2001xc}.
This brings an unknown error into the game --- the error presented in \cite{BG} is
an estimate by the authors obtained by taking the uncertainty of the amplitudes
equal to one half of the two-loop contributions.

Both the Dalitz plot parameter discrepancy and this drawback of ChPT
affecting the predictivity of the chiral computation contributed to the development of
alternative approaches, among others the dispersive methods \cite{KWW,Anisovich:1996tx,Colangelo:2009db}
and the non-relativistic effective field theory (NREFT)
\cite{Colangelo:2006va,Gasser:2011ju,Bissegger:2007yq,Gullstrom:2008sy,Schneider:2010hs}.
In order to understand their relative advantages and disadvantages, let us recall a few basic properties they share. All these approaches are constructed as effective field theories that on the basis of some assumptions (usually represented by some expansion of the amplitude) divide the phase-space of each amplitude into the ``low-energy part'' that is included in the computation and the ``high-energy part'' that is not known or at least less known. At tree-level one simply uses
the amplitudes only in the low-energy region and is not concerned by what lies above
the cutoff. In order to work consistently one needs to introduce a mechanism
that picks up the contributions that contribute with the same importance to a
given order, usually represented by a power-counting. Then,
when computing the amplitude to the higher order, one needs to include also loop
contributions (either by means of taking into account loop Feynman diagrams, as
a unitarity contribution, or by any other method), where one has to integrate
also over the high-energy part of the intermediate amplitudes (over higher
momenta of the intermediate virtual particles). By using the ``power-counting mechanism'' or
by adding some further assumptions, part of these contributions are
considered negligible, but there always remains a part that is finite and
unknown and has to be parametrized somehow --- usually there occur new effective
parameters in the model and the old ones are renormalized or shifted. Note that in ChPT
that represents a Lagrangian effective field theory the power-counting mechanism is given by the chiral counting, which also monitors the number of LECs (effectively containing the contribution of the physics above the chiral cutoff
--- the hadronic scale) appearing at a given order.

The importance of the one-loop (and in the recent years, also of the two-loop)
$\pi\pi$ rescattering corrections has led \cite{KWW,Anisovich:1996tx,Colangelo:2009db} to abandon,
in a certain sense, strict chiral counting, instead attempting to obtain the amplitudes with two-pion rescattering
effects formally included to all orders. These approaches employ a restricted version of
unitarity (taking into account just the two-pion intermediate state), in the context of
dispersion relations, the aim being to find a numerical fixed point
solution of them. The mechanism assigning the importance to a given contribution
is therefore based on the assumption that the two-pion rescattering effects are
dominant. In the low-energy part of the amplitudes, the unitarity contribution of the physics above the threshold, where
further intermediate states contribute and where the $S$ and $P$ partial waves
of the considered amplitudes cease to be the dominant ones, are taken into account through
subtraction constants. However, in order to restrict their number to a reasonable amount, one needs to impose some
assumptions on the high-energy region (of both the physical amplitudes and of
the amplitude constructed iteratively by the numerical method). In
\cite{Anisovich:1996tx,Colangelo:2009db} these assumptions are specified by the requirement to have only
four\footnote{\label{ftnt o poctu subtrakci}Here we
classify the number of parameters appearing in the case we take masses of the
charged and the neutral pions equal, $m_{\pi^\pm}=m_{\pi^0}$. Note that from
Sec.~\ref{sec:conwithchpt} it is obvious that in two-loop ChPT results there
occur at least six independent combinations of LECs.} of them.

The methods based on the modified non-relativistic effective field theory
(NREFT)
\cite{Colangelo:2006va,Gasser:2011ju,Bissegger:2007yq,Gullstrom:2008sy,Schneider:2010hs}
implement instead of the usual chiral expansion a
combined expansion in powers of a formal non-relativistic parameter $\epsilon$ and of a
formal partial-wave $\pi\pi$ scattering-characteristics parameter $a$ (representing scattering lengths and higher threshold shape parameters). The amplitude is then computed to the two-loop order in the NREFT Lagrangian
formalism.
The power-counting scheme is therefore based on the non-relativistic expansion
together with the loop expansion (equivalent in this case to the expansion in
the pion scattering parameters). In \cite{Schneider:2010hs} the results are presented
including the orders up to $O(\epsilon^4), O(a\epsilon^5), O(a^2\epsilon^4)$, and
partially also $O(a^2\epsilon^6)$ and $O(a^2\epsilon^8)$. By assuming that the
included orders are dominant, the contribution of the intermediate states other
than the two-pion ones have to be included through four\footref{ftnt o poctu
subtrakci} parameters coming from local interaction terms.

Naturally, the reasonable question that has to be addressed in the future is
whether each set of assumptions (either of ChPT, of the dispersive approaches,
or of the modified NREFT) adequately describes the physics, and whether a possible
drawback in this respect in any of them is paid off by the other advantages it
possesses. The advantages and the disadvantages of these approaches were
nicely summarized in \cite{Kubis-talk}. We emphasize only that NREFT provides
analytic results that are easy to extent beyond the $m_{\pi^\pm}=m_{\pi^0}$ limit, while the
dispersive methods proceed numerically and their extension to full isospin
breaking was never studied. On the other hand, whereas the NREFT expansion in
$\epsilon$ is safe only inside the Dalitz region, the results of the dispersive
approaches should work also in some larger regions beyond it.
Both of them have in common that, in contrast to ChPT, they directly use the
physically measured $\pi\pi$ scattering parameters as inputs, but there remain
four\footref{ftnt o poctu subtrakci} free parameters that have to be fixed
either from matching to ChPT or to experimental data. Moreover, these decays
depend on $R$ or $Q$ merely just through the normalization, which is
factored out in both of methods. Thus, even if those
representations are fitted to experimental data, the determination of $R$ or $Q$ would
still require to match with ChPT at least at one point, thereby fixing the
normalization.

The matching is not an easy task in this context. In addition to the differences in the structures
of these results, since we are matching two different approaches with different
power-counting schemes and assumptions, we need to find the region (or as
discussed above at least one point) and the appropriate orders in both approaches in which
their results are compatible. Nevertheless, thanks to the easy
form of the one-loop ChPT amplitude, and to the fact that the physical regions of
$\eta\to3\pi$ decays are quite small, in both approaches the matching to
one-loop ChPT was obtained (cf.~\cite{Anisovich:1996tx,Colangelo:2009db} and \cite{Schneider:2010hs}).

In conclusion, in order to determine the correct value of $R$ from the $\eta\to3\pi$ decays, one cannot
avoid discussing either the accuracy of the ChPT result for the amplitude and its possible corrections (by correcting the values of the $O(p^6)$ LECs $C_i$ or by inclusion of some higher-order corrections into the ChPT calculation), or the existence of at least one point (or some region) where the current chiral result reproduces well the complete physical amplitude. For instance, the discussion of the influence of the $C_i$s on the results can be studied using directly the ChPT amplitude, but its complexity and its extreme length together with the fact that it includes
many such $C_i$s complicates such an analysis.

In \cite{kknzprep} we have worked out a method using the dispersive relations
and perturbative unitarity (i.e.~a dispersive approach) for the construction
of a representation of the $\eta\to3\pi$ decay amplitudes.
We have given up on including the $\pi\pi$ rescattering contributions to all orders, but
have instead required to obtain an analytic representation and paid care to the assumptions we
are using in the construction, thus ensuring that the ChPT result can be
obtained as a special case of our result. The method is based on very general
principles, unitarity, analyticity, crossing symmetry, and relativistic
invariance, combined with chiral counting. The fact that we require a
representation valid to two loops in the chiral counting picks up the
contributions that have to be included into the computation and tells us that at the low-energy region up to this chiral order all the other effects are taken into account effectively in terms of six subtraction
parameters\footref{ftnt o poctu subtrakci}.

The full strength (and our original motivation) of this method arises when we want to include the isospin breaking induced by the mass differences between mesons belonging to the same isomultiplet (cf.~\cite{kknzprep}). However, even in the case where we consider the leading order in the isospin breaking, for which the two-loop ChPT result is available, the representation constructed by this method can be useful. Thanks to its simple and compact analytic form, to its capabilities to include all the chiral $O(p^6)$ effects important in the kinematic decay region of $\eta\to3\pi$ into those six subtraction parameters, and to its easy correspondence to the ChPT, this representation is helpful when one is addressing the questions we have premised above, namely, whether one can obtain a reasonable agreement in the determination of the Dalitz parameters from experiment and from the NNLO ChPT amplitude with the corrected set of the $C_i$s; how such a change would influence the determination of $R$; possibly also whether there exists another simple way how to solve that disagreement.

In addition, we do not need to work only in such a close connection to the two-loop ChPT amplitude. Our representation is more general than the two-loop ChPT amplitude (simply stated in the way that the values of our parameters need not to be held at the values stemming from the ChPT), based only on the specific chiral orders of the partial waves of the amplitudes (cf.~e.g.\ \cite{nas-proceeding2}). In order to respect such chiral power-counting, we need to distinguish between various orders of our parameters. By weakening this requirement and by a simple change of their interpretation we
can perform a partial resummation that mimics a part of the previous dispersive approaches. By that we have therefore replaced the assumptions represented by the chiral counting with the assumption that the contributions we have included by this resummation are the dominant one.

In any case, such representation is suitable for fitting the experimental data. We can thus change completely the strategy and instead of trying to correct the amplitude stemming from ChPT, we use our representation as a parametrization of the data, from which we can compute the value of $R$. However, as was discussed above also in this case, we need to fix the normalization from ChPT. For that end we need to find a region where the chiral expansion of the amplitude converges fast, where the two-loop ChPT amplitude reproduces the physics well. Thanks to the form of our representation and its simple
connection to ChPT the analytic dispersive method is helpful also in this analysis, resulting with the recipe for such matching.

We want to stress that in \cite{kknzprep} the inclusion of the isospin-breaking corrections stemming from $m_{\pi^\pm}\neqm_{\pi^0}$ is presented, but as we discuss in Sec.~\ref{sec:neutral} the current experimental data do not yet allow to
perform the isospin-breaking analysis. Thanks to the planned improvement in the neutral decay measurements
(cf.~\cite{Unverzagt:2009vm}) we should add that this possibility is just behind
the
corner. In the paper we therefore work in the limit $m_{\pi^\pm}=m_{\pi^0}$ with the
exception of a few discussions of the effects appearing beyond this limit. This
discussion in full detail is however planned in our next paper \cite{kknzprep2}.

Through relation (\ref{M0Mxiso}) this limit connects the charged
$\eta\to\pi^+\pi^-\pi^0$ decay with the neutral one, $\eta\to3\pi^0$. Using our
representation on the basis of the above mentioned analyses of the charged KLOE
data \cite{KLOE-etac} (the most precise measurement of this process that exists), we can therefore determine the values of the neutral Dalitz decay parameter $\alpha$ and discuss its connection to the direct
neutral measurements (from Table~\ref{tabulka Dalitz eta0}).

The plan of our paper is as follows. After recalling some basic properties of
the amplitudes of these decays and introducing our notation in
Sec.~\ref{sec:notation}, we recall in Sec.~\ref{sec:chpt} the ChPT computation
of $\eta\to3\pi$ with the special emphasis on the contribution of $O(p^6)$ LECs
to the Dalitz plot parameters. From that analysis there follow a few combinations
of observables that should be (at least in the first approximation) safe from the
incorrect determination of these LECs. In Sec.~\ref{sec:dispersive construction}
we present the dispersive construction of our representation for the $\eta\to3\pi$ decay amplitudes. Section~\ref{sec:conwithchpt} discusses the connection between our
representation and the ChPT result. Section~\ref{sec:charged} is then devoted to
the numerical analysis of the charged decay. We start with the determination of the values of our parametrization for NNLO ChPT. Then, inspired by the result of Sec.~\ref{sec:chpt}, we study the influence of changing the $C_i$s in the NNLO ChPT amplitude in order to reproduce the charged KLOE data \cite{KLOE-etac} on the values of the physical observables we are interested in. Then in Sec.~\ref{sec:res} we perform the analysis in which the values of the dispersive parameters are set by KLOE and only the normalization is
determined from ChPT. In that section we present also the procedure of the matching that should reduce the uncertainty coming from the chiral expansion of the amplitudes. In Sec.~\ref{sec:neutral} we use these analyses also for the determination of the neutral Dalitz decay parameter $\alpha$ and discuss briefly the determination of the ratio of the neutral and the charged decay width.
Finally, our conclusions can be found in Sec.~\ref{sec:conclusions}.
We have devoted two appendices to the discussion of properties of the kinematic
functions appearing in our dispersive representation.

\section{Basic properties}\label{sec:notation}

\subsection{Kinematics and notation}
We are interested in two decay channels of $\eta$, the charged one
$\eta\to\pi^+\pi^-\pi^0$ and the neutral one $\eta\to3\pi^0$, generically denoted as
\begin{equation}
 \eta(P) \to \pi^a (p_1) \pi^b (p_2) \pi^0 (p_3)\,.
\end{equation}
The amplitudes of these decay processes can be obtained by analytic continuation of the
amplitudes of the corresponding scattering process
\begin{equation}
 \eta(P) \pi^0(p) \to \pi^a(p_1) \pi^b(p_2)
\end{equation}
by taking $p=-p_3$.
In the following sections, we use the usual Mandelstam variables. In the
scattering region they are defined by
\begin{equation}
s=(P+p)^2,\qquad t=(P-p_1)^2,\qquad u=(P-p_2)^2 ,
\end{equation}
while in the decay region we take
\begin{equation}
s_j=(P-p_j)^2.
\end{equation}
These variables satisfy
\begin{equation}\label{soucet Mandelstamu decay}
s+t+u=3s^c=s_1+s_2+s_3,
\end{equation}
where $s=t=u=s^c$, with
\begin{equation}
3s^c=m_\eta^2 + 2 m_{1,2}^2 + m_3^2
\label{cdp}
\end{equation}
 corresponds the center of the Dalitz plot.
Here $m_3 = m_{\pi^0}$, whereas $m_{1,2}=m_{\pi^0}$ for $\eta\to3\pi^0$, while $m_{1,2}=m_{\pi^\pm}$ for $\eta\to \pi^+ \pi^- \pi^0$.
Up to a convention-dependent phase factor, the crossing relation then means a substitution of the
variables $(s,t,u)$ by $(s_3,s_1,s_2)$, together with the appropriate analytic
continuation from the scattering to the decay region. Bearing this in mind, we can therefore
interchange freely between these two sets of variables.

The constraints (\ref{soucet Mandelstamu decay}) tell us that just two of the kinematic variables are independent.
We can choose them to be, for instance, $s_3=s$
and $s_1=t$. The plot of the dependence of the decay amplitudes on these variables is called Dalitz plot. The physically
allowed kinematical regions for the different crossed amplitudes are constrained by
kinematical limits arising from the condition that the energy of a real particle
has to be at least equal to its rest mass.
Therefore, for a decay process the variable $s_3$ is bounded by
\begin{equation}
(m_1+m_2)^2\le (p_1+p_2)^2=s_3=(k-p_3)^2\le (m_\eta-m_3)^2 ,
\end{equation}
whereas for a scattering in the $s$-channel
\begin{equation}
s\ge(m_\eta+m_3)^2.
\end{equation}
For a fixed value of $s_3$, we obtain bounds for the physical values of $s_1$ (and similarly for $s_2$),
$s_1^{-}(s_3)\le s_1\le s_1^{+}(s_3)$ with
\begin{equation}
s_1^\pm(s_3)=\frac{3s^c-s_3}{2}-\frac{\Delta_{\eta3}\Delta_{12}}{2s_3}\pm\frac{\lambda_{\eta3}^{1/2}(s_3)
\lambda_{12}^{1/2}(s_3)}{2s_3} \,,
\end{equation}
where ($i,j = \eta,1,2,3$)
\begin{align}
\Delta_{ij} &=m_{i}^{2}-m_{j}^{2},\\
\lambda_{ij}(s_3)&=\left(s_3-(m_{i}+m_{j})^2\right)\left(s_3-(m_{i}-m_{j})^2\right). \label{lambda}
\end{align}
Since for both cases under study $m_1=m_2$, the bounds simplify to
\begin{equation}
s_1^\pm=\frac12\left(3s^c-s_3\pm\lambda_{\eta3}^{1/2}(s_3)
\sigma(s_3)\right) ,
\end{equation}
with
\begin{equation}
\sigma(s_3)=\sqrt{1-\frac{4m_{1,2}^2}{s_3}}\,. \label{sigma}
\end{equation}

As was recalled in the Introduction, the amplitudes of the processes $\eta \to 3\pi$, $\mathcal{A}(s,t,u)$, are proportional to the difference of $m_u$
and $m_d$. We therefore pull out this factor, defining
\begin{equation}\label{ARM}
\mathcal{A}(s,t,u) = \frac{\sqrt3}{4}\frac{1}{R}\frac{1}{F_\pi^2} \mathcal{M}(s,t,u)\,,
\end{equation}
where the ratio $R$, which is defined in~(\ref{RQ}), measures the relative violations of $SU(3)$ and of $SU(2)$, and $F_\pi$ is the physical pion decay rate (in our numerical analyses of Sec.~\ref{sec:charged} in order to be in correspondence with \cite{BG}, we take\footnote{Recent analyses (e.g.\ \cite{Kampf:2009tk}) indicate a slightly smaller physical value, $F_\pi=(92.22\pm0.07)\,\mathrm{MeV}$. In order to fully include this change into our computation, redoing of the analysis \cite{BG} with new values for $F_\pi$ and for the pseudoscalar masses would be necessary. Note, however, that a mere change of this value just in this definition leads to a shift of the value of $R$ of about 0.4\%, which is negligible with respect to the other sources of errors occurring in the presented results.} $F_\pi=92.4\,\mathrm{MeV}$).
In accordance to the notation introduced in our general
paper \cite{kknzprep}, when the distinction becomes necessary, the quantities associated to the charged ($\eta\to \pi^+\pi^-\pi^0$) or neutral decay
($\eta\to3\pi^0$) are denoted with the subscript $x$ or $0$, respectively.

In this paper we work mainly to lowest order in the isospin breaking, i.e. we consider the case where all isospin breaking is contained already
in the normalization prefactor $R$ from \eqref{ARM}, and the rest of the amplitude is computed in the isospin limit.
Then due to the isospin structure, the amplitudes $\mathcal{M}(s,t,u)$ are related by
\begin{equation}\label{M0Mxiso}
{\cal M}_0 (s,t,u) = -{\cal M}_x(s,t,u) - {\cal M}_x(t,u,s) - {\cal M}_x(u,s,t)
\end{equation}
(the minus sign is due to the Condon and Shortley phase convention) and in both $\mathcal{M}_0(s,t,u)$ and
$\mathcal{M}_x(s,t,u)$ there appears only one pion mass $m_\pi$. It is why we refer to this case as the $m_{\pi^\pm}=m_{\pi^0}$ limit,
 or more loosely as the isospin limit. However, when making comparisons with the ChPT calculation of \cite{BG}, we use exactly the same values for $\eta$ and $\pi$ masses as were used there\footnote{Note that the value used $m_\eta = 547.3~\text{MeV}$ differs slightly from the current PDG value \cite{PDG}.}, and take for the isospin mass $m_\pi$ in each process a different value --- in the case of the charged $\eta$ decay we take $3 m_\pi^2 =2 m_{\pi^\pm}^2 + m_{\pi^0}^2$, whereas $m_\pi=m_{\pi^0}$ in the case of the neutral decay. So defined $m_\pi$ has the advantage that in both processes we reproduce exactly the physical location of the center of the Dalitz plot and reproduce almost exactly the physical value of the normalization $Q_\eta$ of the Dalitz variables $x$ and $y$ from (\ref{kinematic variables eta})
below. When computing the integrations over the phase space used for setting the normalization from the measured decay rate, we employ again the physical $m_{\pi^\pm}$ and $m_{\pi^0}$ masses for the determination of the phase space.

\subsection{Dalitz plot parametrization}\label{Dalitz plot parametrization}
The standard parametrization of a decay process $P\to3\pi$ is called a
Dalitz plot parametrization (cf.~\cite{PDG}). It is a polynomial
expansion of $|\mathcal{A}(s_1,s_2,s_3)|^2$ around the center of the Dalitz plot.
The parameters are usually normalized in order to be
dimensionless.
The variables of standard use for the charged $\eta$ decay are then
\begin{equation}\begin{split}
x=\sqrt{3}\,\frac{T_1-T_2}{Q_\eta}=\frac{\sqrt{3}\,(s_2-s_1)}{2m_\eta
Q_\eta}\,,\qquad
y=\frac{3T_3}{Q_\eta}-1=\frac{3}{2m_\eta
Q_\eta}\left((m_\eta-m_{\pi^0})^2-s_3\right)-1,
\end{split}\label{kinematic variables eta}
\end{equation}
where $T_j$ is the kinetic energy of the $j$-th pion in the $\eta$ rest-frame. For the charged decay the energy of the reaction $Q_\eta=m_\eta-2m_{\pi^\pm}-m_{\pi^0}$ whereas for the neutral one $Q_\eta=m_\eta-3m_{\pi^0}$. In the case we use in this definition the physical values of the masses, for the charged decay the point $x=y=0$, around which we expand the amplitude, does not coincide\footnote{The $x=y=0$ point is slightly shifted in the $s_3$ direction to $s_3=s^c_x+\frac23(m_{\pi^\pm}-m_{\pi^0})(2m_\eta-m_{\pi^\pm}-m_{\pi^0})$.} exactly with the center of the
Dalitz plot. However, in the isospin limit,
\begin{equation}
y=\frac{3}{2m_\eta Q_\eta}(s^c_x-s_3)
\end{equation}
and the center of the expansion $x=y=0$ matches the center of the Dalitz plot.

The parameters relevant to the decay $\eta\to\pi^+\pi^-\pi^0$ are usually labeled according to
\begin{equation}
\label{dalitz etac}
|{\cal M}_x(s_1,s_2,s_3)|^2=|A|^2\big(1+a y+b y^2+c x+ d
x^2 + e x y+ f y^3+g x^2 y\dots\big)\,,
\end{equation}
where $A$ is the value of the amplitude ${\cal M}_x$ at the point $x=y=0$. Charge conjugation forbids the appearance of terms containing odd powers of $x$ in this expansion, and so $c=e=0$.

\begin{table*}
\begin{center}
\begin{tabular}{c@{\hspace{3.5em}}r@{$\:\pm\:$}l@{\hspace{3em}}r@{$\:\pm\:$}l@{\hspace{3em}}r@{$\:\pm\:$}l@{\hspace{3em}}r@{$\:\pm\:$}l@{\hspace{0.5em}}}
\hline\hline
&\multicolumn{2}{c@{\hspace{3.5em}}}{$a$}&\multicolumn{2}{c@{\hspace{3em}}}{$b$}&\multicolumn{2}{c@{\hspace{3em}}}{$d$}
&\multicolumn{2}{c@{\hspace{0.5em}}}{$f$}\\
\hline
Gormley et al.\cite{Gormley:1970qz} & -1.17\phantom{0} &0.02 & 0.21\phantom{0} &0.03 & 0.06\phantom{0} &0.04
&\multicolumn{2}{c@{\hspace{0.5em}}}{}\\
Layter et al.\cite{Layter:1973ti} & -1.08\phantom{0} &0.014 & 0.034&0.027 &
0.046&0.031&\multicolumn{2}{c@{\hspace{0.5em}}}{}\\
Crystal Barrel \cite{CrystalBarrel-etac} & -1.22\phantom{0} &0.07 & 0.22\phantom{0} &0.11 &
0.06\phantom{0} &0.04&\multicolumn{2}{c@{\hspace{0.5em}}}{}\\
KLOE \cite{KLOE-etac} & -1.090&0.020 & 0.124&0.012 & 0.057&0.017 & 0.14\phantom{0} &0.02\\
\hline
ChPT NNLO \cite{BG} & -1.271&0.075 & 0.394&0.102 & 0.055&0.057 &
0.025&0.160\\
\hline\hline
\end{tabular}
\caption{Values of the Dalitz plot parameters of the $\eta\to\pi^+\pi^-\pi^0$ decay
coming from various experimental and theoretical
determinations.}\label{tabulka Dalitz etac}
\end{center}
\end{table*}

The values of the parameters obtained by various experiments are listed in Table~\ref{tabulka Dalitz etac}. These values are compared with the NNLO calculation in ChPT \cite{BG}. All of
the experiments find the values of $c$ and $e$ compatible with zero.
From the table it is obvious that the precision of the determination from KLOE
\cite{KLOE-etac} exceeds significantly the precision of all the other experiments, which
are more than ten years older. It is also up to now the only experiment that has
determined the parameter $f$ with a reasonable precision.

At this point let us also mention the \emph{linear\/} Dalitz parametrization for the amplitude itself (cf. Appendix A of \cite{BG}):
\begin{equation}\label{lindalnab}
 {\cal M}_x(s,t,u) = A ( 1 + \bar a y + \bar b y^2 + \bar d x^2  + \bar f y^3 + \bar g
x^2 y + \ldots)\,,
\end{equation}
where the parameters $\bar{a}, \bar{b}, \dots$ can now be complex in general. (We have already omitted the terms violating the charge conjugation symmetry of the amplitude.) The parameters of (\ref{dalitz etac}) can be expressed in terms of these \emph{linear\/} Dalitz parameters --- the relations are simple to obtain by squaring (\ref{lindalnab}) and by comparing the terms with the same powers of $x$ and $y$.

At leading order, the parametrization of the $\eta\to3\pi^0$ differential decay rate depends only on the
kinematical variable
\begin{equation}
z=\frac{3}{2m_\eta^2(m_\eta-3m_\pi)^2}\sum_j(s_j-s^c_0)^2
=\frac{3}{2m_\eta^2(m_\eta-3m_\pi)^2}
\left(s_1^2+s_2^2+s_3^2-3(s^c_0)^2\right),
\end{equation}
which denotes the distance from the center of the Dalitz plot, normalized to one at
the edge of the decay region. However,
higher orders corrections do not preserve this accidental rotational symmetry, and we
need again $x$ and/or $y$ from relations (\ref{kinematic variables eta}). Note that the relation
\begin{equation}
z=x^2+y^2
\end{equation}
holds. The Dalitz plot parametrization for this process reads
\begin{equation}\label{dalitz eta0}
\frac{|{\cal M}_0(s_1,s_2,s_3)|^2}{|{\cal M}_0(s^c_0)|^2}=1+2\alpha z+2\beta
y(3z-4y^2)+\gamma z^2+\dots\,.
\end{equation}
The factor of $2$ in front of $\alpha$ and $\beta$ is a mere convention to stress the
connection with the direct \emph{linear\/} Dalitz parametrization of the amplitude
itself (see (\ref{lindalnab}) above). For a better visualization of the violation of the rotation symmetry in $(x,y)$- plane at higher orders, it is convenient to introduce the polar coordinates (cf. also \cite{Schneider:2010hs}), $x=\rho
\cos\phi$, $y=\rho\sin\phi$ with distance $\rho^2=z$, for which we have
$y(3z-4y^2) = \rho^3 \sin(3\phi)$.

\begin{table}
\begin{center}
\begin{tabular}{c@{\hspace{3em}}r@{$\:\pm\:$}l@{\hspace{1em}}}
\hline\hline
&\multicolumn{2}{c@{\hspace{1em}}}{$\alpha$}\\
\hline
Crystal Barrel \cite{CrystalBarrel-eta0}&-0.052\phantom{0}&0.020\\
Crystal Ball \cite{CrystalBall}&-0.031\phantom{0}&0.004\\
WASA/CELSIUS \cite{Celsius}&-0.026\phantom{0}&0.014\\
WASA/COSY \cite{WASAatCOSY}&-0.027\phantom{0}&0.009\\
Crystal Ball {@} MAMI-B \cite{MAMI-B} &-0.032\phantom{0}&0.003\\
Crystal Ball {@} MAMI-C \cite{MAMI-C} &-0.0322&0.0025\\
KLOE \cite{KLOE-eta0} &-0.0301&0.0050\\
\hline
ChPT NNLO \cite{BG} &0.013\phantom{0}&0.032 \\
\hline\hline
\end{tabular}
\caption{Experimental and theoretical values of the slope parameter $\alpha$ of the
$\eta\to3\pi^0$ decay.}\label{tabulka Dalitz eta0}
\end{center}
\end{table}

Various experimental and theoretical determinations of the parameter
$\alpha$ are given in Table~\ref{tabulka Dalitz eta0}. Note the sign discrepancy
between the ChPT determination (with however large error bars) and the experimental measurements, which we
will briefly address in Sec.~\ref{sec:neutral}. Up
to now, no experiment has so far published any constraint on the other parameters, such as
$\beta$.

In the case we work to first order in isospin breaking, the isospin relation (\ref{M0Mxiso}) translates into the following
relations between the neutral Dalitz parameters and the parameters of the \emph{linear\/} parametrization (\ref{lindalnab}) (cf.\ again Appendix~A of \cite{BG})
\begin{align}
\alpha&= \frac12 (\re \bar b + \re \bar d),\label{iso-alpha-bar}\\
\beta&= \frac14 (\re \bar g - \re \bar f).\label{iso-beta-bar}
\end{align}
They can be rewritten in terms of Dalitz parameters of the charged decay. However, there still remains a dependence on the imaginary parts of the \emph{linear\/} parameters,
\begin{align}
\alpha&=\frac14\left(b+d-\frac{a^2}{4}\right)-\frac{\left(\im \bar{a}\right)^2}{4},\label{iso-alpha}\\
\beta&=\frac18\left(g-f+\frac{a}{2}(b-d)-\frac{a^3}{8}\right)+\frac{\im \bar{a}}{4}
\left(\im \bar{b}-\im \bar{d}-a \frac{\im\bar a}{4}\right).
\label{iso-beta}
\end{align}

\subsection{Adler zero}\label{sec:adler}
The isospin-breaking part of the QCD Hamiltonian density (\ref{eta hamiltonian QCD}) can be written as (in this subsection $\lambda^a$ are Gell-Mann SU(3) matrices)
\begin{equation}
\mathcal{H}^{IB}_{\mathrm{QCD}}(x)=(m_d-m_u)S_3(x),
\end{equation}
where
\begin{equation}
S_3(x)=\frac 12\overline{q}(x)\lambda^3 q(x).
\end{equation}
Therefore, to first order in $m_d-m_u$, the amplitudes of the isospin-breaking processes that are described by this Hamiltonian
are connected to Green functions with one insertion of zero momentum scalar density $S_3$ (calculated in the limit $m_u=m_d$). In
the $SU(2)\times SU(2)$ chiral limit $m_u=m_d=0$, pions are genuine Goldstone bosons. For the corresponding amplitudes $\langle f+\pi
^a(p)|S_3(0)|i\rangle $ with a pion in the final state, we can thus derive the $SU(2)\times SU(2)$ soft-pion theorem in the general form
\begin{equation}
\langle f+\pi ^a(p)|S_3(0)|i\rangle \rightarrow 0\quad\mathrm{for}\quad
p\rightarrow 0.
\end{equation}
The derivation of the theorem proceeds in the usual way, except that now, because of the insertion
of $S_3(0)$ transforming under the axial $SU(2)$ rotation nontrivially as
\begin{equation}
\delta_5^aS^3(0)=\frac{\mathrm{i}}2\overline{q}(x)\{\lambda^3,\lambda
^a\}\gamma_5 q(x),
\end{equation}
it only  holds provided $a=1,2$. (For $a=3$ there occurs an additional contribution from $\delta_5^3 S^3(0)$, which does not vanish.)

For the decay $\eta \to \pi^{+}(p_{+})\pi^{-}(p_{-})\pi^0(p_0)$, this means that the amplitude $\mathcal{M}_x(s,t,u)$ defined in (\ref{ARM}) vanishes in the $SU(2)\times SU(2)$ chiral limit for either $p_{+}=0$ or $p_{-}=0$, i.e.\ it develops two Adler zeroes \cite{Adler:1965ab,Adler:1965ga} $s=u=0$, $t=m_\eta^2$ and $s=t=0$, $u=m_\eta^2$.
As a consequence, expanding the amplitude $\mathcal{M}_x(s,t,u)$ beyond the $SU(2)\times SU(2)$ chiral limit in the independent variables $s$ and $\left(t-u\right) $ around the points where $s=0$, $(t-u)^2=m_\eta^4$, or more specifically around the points
\begin{equation}\begin{split}
\overline{s} &=0, \\
\overline{t} &=\frac 12(m_\eta ^2+3m_\pi ^2\pm m_\eta ^2), \\
\overline{u} &=\frac 12(m_\eta ^2+3m_\pi ^2\mp m_\eta ^2),
\end{split}\end{equation}
according to (here we use the $t\leftrightarrow u$ symmetry of the amplitude)
\begin{equation}
\mathcal{M}_x(s,t,u)=\sum_{i,j\geq 0}c_{ij}s^i((t-u)^2-m_\eta ^4)^j,
\label{A_zero_expansion}
\end{equation}
we can restate the above $SU(2)\times SU(2)$ theorem in the form
\begin{equation}
\lim_{m_u,m_d\rightarrow 0}c_{00}=0.
\end{equation}
Since the position of the Adler zero is determined up to $O(m_\pi^2)$ corrections, an analogous statement remains true also for similar expansion coefficient $c_{00}^{(\xi,\zeta)}$ corresponding to an expansion around the
points with $s=\xi m_\pi^2$, $(t-u)^2=\left( m_\eta ^2+\zeta m_\pi ^2\right)^2$, namely,
around the points
\begin{equation}\label{adlerxi}\begin{split}
s_{(\xi,\zeta)} &=\overline{s}+\xi m_\pi ^2, \\
t_{(\xi,\zeta)} &=\overline{t}\pm \frac 12(\zeta\mp \xi)m_\pi ^2, \\
u_{(\xi,\zeta)} &=\overline{u}\mp \frac 12(\zeta\pm \xi)m_\pi ^2,
\end{split}\end{equation}
where $\xi,\zeta$ are reasonably small and behave as $O(1)$ for $m_u,m_d\rightarrow 0$.
For the value of the amplitude at these points we therefore obtain
\begin{equation}
\mathcal{M}_x(s_{(\xi,\zeta)},t_{(\xi,\zeta)},u_{(\xi,\zeta)})=c_{00}^{(\xi,\zeta)}=O(m_\pi^2)
\end{equation}
and its absolute value is expected to be numerically small.

Note, however, that the remaining coefficients $c_{ij}$ are not protected by such a factor $m_\pi ^2$, and the same is true
also for the value of the amplitude at points far from $(\overline{s},\overline{t},\overline{u})$, where
$\mathcal{M}_x(s,t,u)$ can be enhanced by a factor $m_\eta^2/m_\pi^2$ with respect to
$\mathcal{M}_x(\overline{s},\overline{t},\overline{u})$.
Note also that a small numerical value of $c_{00}$ (or $c_{00}^{(\xi,\zeta)}$ in general) does not
necessarily imply that its chiral expansion shows better convergence than
the one of any other $c_{ij}$, in the sense that for the ratio of two subsequent
corrections the relation
\begin{equation}
\frac{c_{00}^{(n+1)}}{c_{00}^{(n)}}\ll \frac{c_{ij}^{(n+1)}}{c_{ij}^{(n)}}
\end{equation}
does not necessarily hold.

\subsection{Isospin violation and cusp}\label{sec:cusp}
In the case we go beyond the first order in the isospin breaking, in addition to the complications that the two $\eta\to3\pi$ decay amplitudes are no longer connected by (\ref{M0Mxiso}), and that the expressions are more complicated due to the fact that there appear two different masses of pions, in the processes with two neutral pions in the final state there occurs an interesting phenomenon called cusp. This effect is caused by different charged and neutral pion masses and is connected with the contributions of $\pi^+\pi^-$ intermediate states rescattering back to $\pi^0\pi^0$. Such  a state generates a square root singularity, which resides at $4m_{\pi^\pm}^2$, lying above the physical threshold, $4m_{\pi^0}^2$, and the unitarity cusp is a result of the interference between the part of the amplitude containing this singularity and the
rest without it.

It is obvious that the cusp emerges only in the case when isospin breaking is included also in ${\cal M}(s,t,u)$ and that its
strength is sensitive to $\pi^+\pi^-\to\pi^0\pi^0$ scattering at threshold (mainly to the scattering length of this process). This property can be used for a determination of the scattering lengths from the measurement of the cusp \cite{Cabibbo:2004gq,Batley:2005ax,Abouzaid:2008js}.

Let us try to estimate the relative sizes of the cusps in various processes where a pseudoscalar, namely $K^+$, $K_L$ or $\eta$, decays into three pions. (This discussion is inspired by \cite{CI} and \cite{DiLella-talk}.) Because the pion rescattering part will be approximately the same for all the processes, we may consider the notion of ``visibility'' of the cusp in these processes by comparing the relative ratios between the cusps and the regular parts of the amplitudes,
\begin{equation}
\gamma(P^c)=\left.\frac{\kappa_c |\mathcal{A}_{P^c\to\pi^+\pi^-\pi^c}|\
|\mathcal{A}_{P^c\to\pi^0\pi^0\pi^c}|}{|\mathcal{A}_{P^c\to\pi^0\pi^0\pi^c}|^2}
\right\rvert_{\text{cusp}}
=\left.\kappa_c\frac{|\mathcal{A}_{P^c\to\pi^+\pi^-\pi^c}|}{|\mathcal{A}_{P^c\to\pi^0\pi^0\pi^c}|}
\right\rvert_{\text{cusp}},
\end{equation}
where $|\mathcal{A}|$ is the absolute value of the matrix element of the indicated
process and $\kappa_c$ is a multiplicity factor corresponding to that process,
equal to $2$ in the case the decaying particle is charged (two possible
$\pi^+\pi^-\to\pi^0\pi^0$ scatterings are then possible), and to $1$ in the other
cases. These ratios have to be evaluated at the cusp point $s=4m_{\pi^\pm}^2$.

Using the measured relative decay rates and the values of Dalitz parameters from
\cite{PDG}, we obtain for these processes,
\begin{equation}
\gamma(K^\pm)\sim7.3,\qquad \gamma(K_L)\sim0.45,\qquad \gamma(\eta)\sim0.93.
\end{equation}
From that we can estimate that the effect is approximately 16 (8) times more pronounced in the $K^\pm$
decay with respect to $K_L$ $(\eta)$ decay.

First indications of the cusp effect were already observed also in the $\eta\to3\pi^0$ decay (cf.\ e.g.\ \cite{MAMI-C}). This effect however appears at the edge of the decay region\footnote{In the $(x,y)$ plane, the cusp is located on the segment
$y=y(4m_{\pi^\pm}^2)\sim0.773$ and on two other segments obtained by $s\leftrightarrow t$ and $s\leftrightarrow u$ (i.e.\ obtained by rotation of the original one by $\pm120^\circ$ around the center of the Dalitz plot). Its position thus does not respect the accidental rotation symmetry, and depending on its direction in the $(x,y)$ plane, the corresponding value of $z$ changes from $0.597$ to $0.883$ as $\frac{0.597}{\cos(\phi-\phi_0)}$, with $\phi_0=0^\circ,\pm120^\circ$.} and is therefore not simple to measure.

For the time being, because of this lack of data, we shall not  pursue the discussion about the cusp here (even though our representation describing also this effect is prepared \cite{kknzprep, nas-proceeding1}). Instead we will work in the strict isospin limit beyond the trivial order at which $\eta\to3\pi$ decay is forbidden, i.e. ${\cal M}_x$ or ${\cal M}_0$ is taken in the isospin symmetry limit.

\section{Chiral perturbation theory}\label{sec:chpt}
Let us briefly recapitulate the ChPT calculation of $\eta\to3\pi$ decays.
As was discussed in the Introduction, direct electromagnetic corrections to these processes are very small, and thus they
proceed mainly through the part~(\ref{eta hamiltonian QCD}) of the QCD Lagrangian. The leading order (LO) calculation was performed in \cite{Bell:1996mi,Cronin:1967jq,Osborn:1970nn}, which in our notation\footnote{In this work for the various chiral orders we follow the convention of \cite{BG}, where the amplitudes are at a given order simplified using Gell-Mann-Okubo relations, physical decay constants and physical pseudoscalar masses and the so induced differences are included into higher orders.} takes the very simple form
\begin{equation}
{\cal M}^{(2)}(s,t,u) = \frac{4}{3} m_\pi^2 - s.
\end{equation}
The next-to-leading order (NLO) was provided in~\cite{Gasser:1984pr}. Its form is discussed in Sec.~\ref{sec:etaoneloop} below. The $O(p^6)$ corrections were studied quite recently in~\cite{BG}.
From these three successive orders one can see that $\eta\to3\pi$ thus represents a case where the chiral corrections are large~\cite{Bijnens:2006zp}. Moreover, it seems that also the two-loop ChPT result supplemented with the existing LECs determination of \cite{BG} is not working very well as we have demonstrated on the example of Dalitz parameters in the Introduction.

\subsection{Contribution of the constants \texorpdfstring{$C_i$}{Ci} to Dalitz parameters}
\label{sec:Ci}
In the NNLO result there occurs a great deal of $O(p^6)$ low-energy constants $C_i$ which are only estimated from resonance saturation. Hence, the first question one has to ask is whether the discrepancy with experiment cannot be accounted for by the unsatisfactory knowledge of the $O(p^6)$ low-energy constants.

Let us thus study the contribution of $O(p^6)$ LECs ($C^r_i$) to the Dalitz parameters of the individual decay modes.
There are several possibilities how to determine these parameters from the computed amplitude ${\cal M}(s,t,u)$. For instance, we can expand $|{\cal M}(s,t,u)|^2$ to the order $O(p^6)$, and then make the Taylor expansion at the center of the Dalitz plot. This would result in the linear dependence of the Dalitz parameters on the $C_i$s. Provided we did not chiral expand first and instead made a fit of the modulus squared of the complete $O(p^6)$ amplitude to the Dalitz parametrization (as it was done in \cite{BG}), we would get a more complicated dependence on the $C_i$s including also quadratic and mixed terms.
Such contributions should be, however, suppressed by the chiral counting. Nevertheless, they can bring sizable changes in the final numerical predictions. In order to obtain the linear contribution only, we follow the first procedure.

We start with the neutral decay mode. The explicit dependence of $\alpha$ on the $C_i$s was already given in~\cite{Schneider:2010hs},
\begin{equation}
 \alpha_C = \frac{16 m_\eta^2 (m_\eta-3 m_\pi)^2}{3F_\pi^4}\, {\cal C}_\alpha
\end{equation}
with
\begin{equation}
{\cal C}_\alpha = C^r_5 + C^r_8 + 3 C^r_9 + C^r_{10} - 2 C^r_{12} + 2 C^r_{22}
+ 3 C^r_{24} + C^r_{25}.
\end{equation}
Further, by a careful investigation of the $O(p^6)$ polynomial of the amplitude calculated in \cite{BG},
we realize that there is no contribution of the $C_i$s to the second neutral Dalitz parameter (it is connected with relation (\ref{iso-beta-bar}) and the fact that in the charged decay $\bar f_C = \bar g_C$ as stated below),
\begin{equation}
\beta_C =0.
\end{equation}

In the case of the charged decay we summarize first the contributions of the $C_i$s to the \emph{linear} coefficients $\bar a$, $\bar b$, $\bar d$, $\bar f$, $\bar g$ defined in~(\ref{lindalnab}) that are directly connected with the amplitude. These parameters can be in general complex but since we deal only with the linear contribution of the $C_i$s, they contribute only to their real values.

By a simple algebra one obtains
\begin{equation}
\bar a_C = - \frac{8 m_\eta (m_\eta-3 m_\pi)}{3 F_\pi^4 (m_\eta^2-m_\pi^2)}\,
\left(m_\eta^4 {\cal C}_{a}^\eta -6 m_\eta^2 m_\pi^2 {\cal C}_{a}^{\eta\pi}
-3 m_\pi^4 {\cal C}_{a}^\pi\right),
\end{equation}
where we have slightly more complex structure
\begin{align}
\begin{split}
{\cal C}_{a}^\eta &= C^r_1-2C^r_3-6C^r_4-2C^r_5-3C^r_6+C^r_8+C^r_{10}
+6C^r_{11}-20C^r_{12}-18C^r_{13}+9C^r_{14}-36C^r_{16}+9C^r_{17}
+45C^r_{18}\\
&-81C^r_{19}-54C^r_{20}+8C^r_{22}+9C^r_{24}+C^r_{25}
-18C^r_{27}-36C^r_{28}-54C^r_{31}-54C^r_{32}-108C^r_{33},
\end{split}\\
\begin{split}
{\cal C}_{a}^{\eta\pi}&= C^r_1-2C^r_3-2C^r_4+C^r_5+2C^r_6-C^r_{10}-4C^r_{11}
-4C^r_{12}+12C^r_{13}-4C^r_{14}-3C^r_{15}-12C^r_{16}+2C^r_{17}
+6C^r_{18}\\
&+4C^r_{22}+3C^r_{24}-C^r_{25}-2C^r_{26}-6C^r_{27}
-12C^r_{28}+4C^r_{29},
\end{split}\\
\begin{split}
{\cal C}_{a}^\pi &= C^r_1- 2C^r_3+ 2C^r_4 + 3C^r_6 + 3C^r_8 - 3C^r_{10} - 6C^r_{11}
+ 12C^r_{12} + 18C^r_{13}+ 3C^r_{14} - 6C^r_{15} + 12C^r_{16} - 9C^r_{17}
+ 3C^r_{18}\\
&- 27C^r_{19} - 18C^r_{20} - 3C^r_{24}- 3C^r_{25} - 4C^r_{26}
+ 6C^r_{27} + 12C^r_{28} + 8C^r_{29} - 18C^r_{31} - 18C^r_{32} - 36C^r_{33}.
\end{split}
\end{align}
Similarly, we have
\begin{equation}
\bar b_C = \frac{8 m_\eta^2 (m_\eta-3 m_\pi)^2}{3 F_\pi^4 (m_\eta^2-m_\pi^2)}\,
\left(m_\eta^2 {\cal C}_{b}^\eta + m_\pi^2 {\cal C}_{b}^\pi\right),
\end{equation}
where
\begin{align}
{\cal C}_{b}^\eta &= C^r_1 -2 C^r_3+ 2 C^r_4+4C^r_5+3 C^r_6+C^r_8+6 C^r_9+ C^r_{10}
-6 C^r_{11} -2 C^r_{12}+18 C^r_{13}+2 C^r_{22}+6C^r_{24}+4 C^r_{25},
\\
{\cal C}_{b}^\pi &= 3 C^r_1 -6 C^r_3+6C^r_4-2 C^r_5+3 C^r_6+C^r_8-6 C^r_9-5 C^r_{10}
-6 C^r_{11}+10 C^r_{12}+18 C^r_{13}-10
C^r_{22}-6 C^r_{24}+4 C^r_{25},
\end{align}
\begin{equation}
\bar d_C = -\frac{8 m_\eta^2 (m_\eta-3 m_\pi)^2}{3 F_\pi^4
(m_\eta^2-m_\pi^2)}\,\left(m_\eta^2 {\cal C}_{d}^\eta + m_\pi^2 {\cal C}_{d}^\pi\right),
\end{equation}
with
\begin{align}
{\cal C}_{d}^\eta  &= C^r_1 -2 C^r_3+2 C^r_4+3 C^r_6-3C^r_8-6 C^r_9-3 C^r_{10}
-6 C^r_{11}+6 C^r_{12}+18 C^r_{13}-6 C^r_{22}-6 C^r_{24},
\\
{\cal C}_{d}^\pi  &= 3 C^r_1-6 C^r_3+6 C^r_4+2C^r_5+3 C^r_6+5 C^r_8+6 C^r_9
-C^r_{10}-6 C^r_{11}+2 C^r_{12}+18 C^r_{13}-2 C^r_{22}+6C^r_{24}+8 C^r_{25},
\end{align}
and
\begin{equation}
\bar f_C =\frac{16 m_\eta^3 (m_\eta-3 m_\pi)^3}{3 F_\pi^4 (m_\eta^2-m_\pi^2)}\, {\cal C}_f
\end{equation}
with
\begin{equation}
{\cal C}_f = C^r_1-2 C^r_3 - 2 C^r_4.
\end{equation}

Finally, the contribution of $C_i$s to the parameter $\bar{g}$ is the same as in the case of the parameter $\bar{f}$,
\begin{equation}\label{f rovna se g}
\bar g_C = \bar f_C.
\end{equation}

\newpage

Using these relations in the same spirit as in \cite{Bijnens:2009zd}, we can construct the
combinations of physical (or quasi-physical) quantities which do not depend on any $C_i$:
\begin{enumerate}
 \item $\left(\bar b + \bar d - 2 \alpha \right)\bigr|_C=0$   $\Rightarrow$
    \begin{equation}\label{eqrel1}
      \text{rel}_1 \equiv \left(4 (b+d)-a^2 - 16 \alpha\right)\bigr|_C=0
    \end{equation}
 \item $\left(\bar f - \bar g\right)\bigr|_C= 0$ $\Rightarrow$
    \begin{equation}\label{eqrel2}
     \text{rel}_2 \equiv \left(a^3-4 a b+4 a d+8 f-8 g\right)\bigr|_C =0
    \end{equation}
 \item $\beta\bigr|_C = 0$
 \item ${\cal C}_b^\eta + {\cal C}_b^\pi -{\cal C}_d^\eta - {\cal C}_d^\pi =0$
\end{enumerate}
Let us discuss them in more detail (in a reverse order).
The last expression, of course, does not represent any combination of physical quantities, and so it is not possible to use it directly in connection with any observable. It could be, however, useful for lattice simulations, where one can vary the meson
masses. On the contrary the third relation, stating that the second neutral Dalitz parameter $\beta$ does not depend on any $C_i$, represents a simple possibility, open to future experiments, how to check the ChPT results unaffected by the error stemming from the estimates of $C_i$. Now let us turn our attention to the
relations $\text{rel}_1$ and $\text{rel}_2$. The latter was mentioned in \cite{Bijnens:2009zd}, while the first one was implicitly stated in \cite{BG}.
In fact, $\text{rel}_1$ is a simple consequence of the isospin relation (\ref{iso-alpha}) stating that the $C_i$s do not contribute to $\im\bar{a}$
and should thus be valid for any real contributions to the Dalitz parameters appearing there (not only for the contributions of the $C_i$s).

\begin{table*}
\centering
\begin{tabular}{@{\hspace{1em}}c@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{0.5em}}}
\hline\hline
&\multicolumn{2}{c@{\hspace{2em}}}{KLOE} & \multicolumn{2}{c@{\hspace{2em}}}{ChPT}& \multicolumn{2}{c@{\hspace{2em}}}{ChPT$_g$}& \multicolumn{2}{c@{\hspace{2em}}}{NREFT$_i$}& \multicolumn{2}{c@{\hspace{0.5em}}}{NREFT}\\\hline
$a$
& $-1.09\phantom{0}$  & $0.02$
& $-1.271$  &  $0.075$
& $-1.262$  &  $0.079$
& $-1.279$  &  $0.012$
& $-1.213$  &  $0.014$
\\
$b$
& $0.124$ & $0.012$
& $0.394$  &  $0.102$
& $0.407$  &  $0.108$
& $0.361$  &  $0.021$
& $0.308$  &  $0.023$
\\
$d$
& $0.057$ & $0.017$
& $0.055$  &  $0.057$
& $0.081$  &  $0.089$
& $0.053$  &  $0.003$
& $0.050$  &  $0.003$
\\
$f$
& $0.14\phantom{0}$  & $0.02$
& $0.025$  &  $0.160$
& $0.009$  &  $0.170$
& $0.089$  &  $0.018$
& $0.083$  &  $0.019$
\\
$g$
& \multicolumn{2}{c@{\hspace{1em}}}{$\sim0$}
& \multicolumn{2}{c@{\hspace{1em}}}{$0$}
& $-0.07\phantom{0}$  &  $0.19$
& $-0.043$  &  $0.002$
& $-0.039$  &  $0.002$
\\\hline
$\alpha$
& $-0.030$ & $0.005$
& \multicolumn{4}{c}{$0.013\pm0.016$}
& $-0.024$ &  $0.005$
& $-0.025$ &  $0.005$
\\
$\beta$ &
\multicolumn{2}{c}{}
& \multicolumn{4}{c}{$-0.002\pm 0.025\phantom{\pm}$}
& $-0.004$ & $0.001$
& $-0.004$ & $0.001$
\\\hline\hline
\end{tabular}
\caption{Dalitz parameters for charged and neutral decay mode of $\eta\to 3\pi$
calculated in ChPT at order $O(p^6)$ \cite{BG}, NREFT \cite{Schneider:2010hs} (in addition to their final value given in the ultimate column,
in the penultimate column captioned with NREFT$_i$ we list also the values without isospin-breaking effects included) and as
measured by KLOE \cite{KLOE-etac,KLOE-eta0}. ChPT$_g$ represents a fit redone with $g$
included (but without any higher parameters), cf.\ fit NNLOq in \cite{BG}.}
\label{tab:inputdalitz}
\end{table*}

For the comparison of ChPT results \cite{BG} with the values measured by KLOE \cite{KLOE-etac}
we can use not only the values of the Dalitz parameters summarized in Table~\ref{tab:inputdalitz} but also the combinations of these parameters (\ref{eqrel1}) and (\ref{eqrel2}) that are (at least in a first approximation) $C_i$-independent.
It means that the influence of all physics beyond the pseudoscalar domain (mainly scalar and vector resonances) on these combinations is hidden in $O(p^4)$ LECs $L_i$, which are phenomenologically much better under control than the $C_i$s, thereby providing a clearer theoretical output. We should remark, however, that the independence of all these relations on the $C_i$s occurs only in the case we take $m_{\pi^\pm}=m_{\pi^0}$. Away from this limit the situation can be different and these combinations can still have non-negligible dependence on the $O(p^6)$ LECs.

\begin{table*}
\centering
\begin{tabular}{@{\hspace{1em}}c@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{1em}}}
\hline\hline
&\multicolumn{2}{c@{\hspace{2em}}}{KLOE} & \multicolumn{2}{c@{\hspace{2em}}}{ChPT}& \multicolumn{2}{c@{\hspace{2em}}}{ChPT$_g$}& \multicolumn{2}{c@{\hspace{2em}}}{NREFT$_i$} & \multicolumn{2}{c@{\hspace{1em}}}{NREFT}\\\hline
rel$_1$
& $0.02$  &  $0.12$
& $-0.03$  &  $0.72$
& $ 0.15$  &  $0.79$
& $0.41$  &  $0.12$
& $0.35$  &  $0.13$
\\
$\chi^2$
&\multicolumn{2}{c@{\hspace{2em}}}{} & \multicolumn{2}{c@{\hspace{2em}}}{$0.1$} & \multicolumn{2}{c@{\hspace{2em}}}{$1.2$} & \multicolumn{2}{c@{\hspace{2em}}}{$10$} & \multicolumn{2}{c@{\hspace{1em}}}{$7.6$}
\\\hline
rel$_2$
& $0.12$  &  $0.21$
& $-0.13$  &  $1.4$
& $0.28$  &  $2.1$
& $0.54$  &  $0.19$
& $0.44$  &  $0.20$
\\
$\chi^2$
&\multicolumn{2}{c@{\hspace{2em}}}{} & \multicolumn{2}{c@{\hspace{2em}}}{$1.4$} & \multicolumn{2}{c@{\hspace{2em}}}{$0.6$} & \multicolumn{2}{c@{\hspace{2em}}}{$4.0$} & \multicolumn{2}{c@{\hspace{1em}}}{$2.4$}
\\\hline\hline
\end{tabular}
\caption{Two relations defined in~(\ref{eqrel1}) and (\ref{eqrel2}) compared
using different models. The quantity $\chi^2$ defined in~(\ref{chi2}) was computed for the theoretical frameworks by taking for the experiment the values from $\text{exp}\equiv\text{KLOE}$ as given in the first column.}
\label{tab:howwell}
\end{table*}

The values of these combinations that use the data from Table~\ref{tab:inputdalitz} are presented in Table~\ref{tab:howwell}.
This table indicates that even though the central values of the individual Dalitz parameters determined by ChPT and KLOE
differ, the central values of these two combinations are in a good agreement, which indicates that ChPT is not working at all that badly. Unfortunately the large errors quoted there somehow put down the importance of any conclusions. However, one should bear in mind that the values quoted in Table~\ref{tab:howwell} were computed just using the values and the error bars presented in Table~\ref{tab:inputdalitz} that were attributed mainly from the fitting procedures and are thus strongly correlated. This can affect the positions of the central values by small changes, but primarily the error bars of these combinations are then overestimated. Note that the errors of the Dalitz parameters from ChPT are enhanced also by large systematic uncertainties of the amplitudes entering these fitting procedures. Such uncertainties were caused mainly by uncertainties of the $C_i$s, which should be substantially eliminated in these combinations. We also observe another artifact of the fitting procedure when comparing the values denoted by ChPT and ChPT$_g$ that differ just by the truncation of the Dalitz parametrization at $f$ and $g$, respectively. The combination $\text{rel}_1$ is according to relation (\ref{iso-alpha}) equal to $4(\im\bar{a})^2$, which should be therefore positive. The value denoted by ChPT does not possess this property even though that the value of $4(\im\bar{a})^2$ obtained by a direct fit of the original amplitude in \cite{BG} reproduces well the value given in the column ChPT$_g$.

A similar effect can occur also for the KLOE values since in \cite{KLOE-etac} the value\footnote{Note the different notation of this Dalitz parameter in \cite{KLOE-etac} --- for the parameter denoted in this text by $g$, KLOE uses symbol $h$.} of $g$ was not presented (only its compatibility with zero). As an illustration, we remind the reader that if we added to the values of $a,b,d,f$ measured by KLOE the value of $g=-0.02$ ($g=-0.04$), we would obtain an exact match of the so defined experimental\footnote{Naturally, repeating the KLOE fit with $g$ included would also change the values of the further parameters (cf.\ again the difference between the values from ChPT and ChPT$_g$).} value of $\text{rel}_2$ with the value from ChPT$_g$ (NREFT).

In these two tables we have also studied predictions of NREFT \cite{Schneider:2010hs}. Since that method is built in a different way than ChPT, the combinations of the observables appearing in $\mathrm{rel}_1$ and $\mathrm{rel}_2$ have no special significance there. However, they are still valid combinations of observables and so nothing prevents us from using them for comparison of the predictions from any theory with the experiment. The lesser agreement of NREFT and KLOE in $\mathrm{rel}_1$ was already pointed out in \cite{Schneider:2010hs} in terms of different values of $\im \bar{a}$ stemming from the representation of \cite{Schneider:2010hs} and the one coming from the KLOE measurement and the relation (\ref{iso-alpha}).
Together with the slight inconsistency also in $\text{rel}_2$ depending only on the parameters of the charged decay,
this indicates that there is a problem either on the side of the current determination from the KLOE group or on the side of the NREFT representation.

We conclude this discussion with the statement that a new measurement of the charged Dalitz parameters (possibly taking into account these two relations) would therefore be highly desirable. Before that, we are not able to answer the question whether it is possible to reproduce the physical Dalitz plot distribution with a better determination of the LECs $C_i$ or whether the discrepancy between the ChPT-computed and the experimentally measured distributions has some other origin (slow convergence of the chiral counting, \dots). In addition, should the experimental value confirm the values inconsistent with the predictions of \cite{Schneider:2010hs}, even if one accepts the explanation for the discrepancy of the neutral parameter $\alpha$ proposed in \cite{Schneider:2010hs}, the issue of the discrepancy for the charged parameter $b$ would remain open.

But for now, inspired just by the quite good consonance of the current KLOE and the ChPT values of the $C_i$-independent relations, we would expect that by finding the right values of the $C_i$s we would reproduce (at least partially) better the physical values of the Dalitz parameters. The natural question can arise now whether it would be possible to find
an elaborate determination of such $C_i$s going beyond the crude estimate of the simple resonance saturation model used
in~\cite{BG}.

Let us start with $\alpha$. Its resonance saturation is simpler as there are no vector resonance
contributions. For the simple scalar resonance model used in~\cite{BG} we obtain (cf.~\cite{Schneider:2010hs})
\begin{equation}
 \alpha|_C = \frac{16 m_\eta^2 (m_\eta-3 m_\pi)^2 c_d c_m}{3 F_\pi^2 M_S^4}
\approx 0.005
\end{equation}
that is positive. However, the minimal chiral symmetry breaking introduced in \cite{Cirigliano:2003yq} changes $c_m$ into $c_m + c_d e_m$, and especially for standard hierarchy ($e_m<0$) one can thus produce a negative contribution to $\alpha$.
Using the same numbers as obtained from the phenomenological study in \cite{Cirigliano:2003yq}, where they distinguish two models, one representing the inverted hierarchy (the model called A) and one representing the standard hierarchy (called B), we obtain
\begin{equation}
\alpha|_{C\,[A]} = 0.002,\; \text{or}\quad
\alpha|_{C\,[B]} = -0.0005\,,
\end{equation}
which lead to the final values $\alpha=0.01$ and $\alpha=0.007$, respectively.

The situation for the charged decay mode looks more complicated. The transition from the amplitude to the Dalitz parametrization leads to many mixing terms and the dependence on the $C_i$s is not linear. Even though, as already mentioned, such higher terms are theoretically suppressed by chiral counting,
in practice they can turn out to be more substantial than anticipated (it is true especially for model $B$). In order
to get more reliable results we perform a full fit to the Dalitz distribution in exact correspondence with \cite{BG}, with the exception that we fit a polynomial of the third order (i.e.\ including $g$), which corresponds to ChPT$_g$ in Table~\ref{tab:inputdalitz}. The vector resonance saturation employed here is based on the model and the phenomenology constraints from~\cite{Kampf:2006bn}. The resulting fits to the Dalitz parametrization are summarized in Table~\ref{tab:dalitz}.
\begin{table}
\centering
\begin{ruledtabular}
\begin{tabular}{crccccc}
model & $|A|^2$ & $a$ & $b$ & $d$ & $f$ & $g$ \\\hline
ChPT$_g$
&534
&-1.26
&0.41
&0.081
&0.009
&-0.072
\\
simple
& 516
& -1.39
& 0.47
& 0.10\phantom{0}
& 0.025
&-0.088
\\
model $A$
&  723
& -1.31
& 0.41
& 0.081
& 0.024
&-0.069
\\
model $B$
&  1835
& -1.19
& 0.33
& 0.052
& 0.020
&-0.040
\end{tabular}
\end{ruledtabular}
\caption{The comparison of the impact of the different models for resonance saturation on
the Dalitz parameters. For vector resonances we have taken the model from~\cite{Kampf:2006bn},
different cases for scalar resonances are: the simple one from \cite{BG}, and the models $A$ and $B$ from \cite{Cirigliano:2003yq}.}
\label{tab:dalitz}
\end{table}
It is clear from this table that, as we have anticipated, the $C_i$s have a bigger effect
than expected from mere chiral counting. They also have an impact on the normalization $|A|^2$, which
in the case of model $B$ is far from being negligible. Let us note at this
point a few things concerning the resonance saturation. It is obvious from
Table~\ref{tab:dalitz} that model $B$ would produce an unrealistic increase
of the amplitude (thereby also of $R$ or $Q$). It does not, however, mean that
this model for scalar resonances is disqualified. Higher resonances, representative of the physics beyond the pseudo-Goldstone bosons, contribute to both $L_i$s and $C_i$s (when talking about NNLO). One cannot just
keep their influence on $C_i$s ignoring their presence in $L_i$s and thus merge
inconsistently two models, i.e.\ in our case the model used in ``fit 10'' of
\cite{Amoros:2001cp} and the models $A$ or $B$. One can always try to be as ``harmless'' as
possible with any extension of the simple resonance saturation and try to
preserve the original values of $L_i$s (as was to some extent possible for the chiral
symmetry breaking construction done in \cite{Kampf:2006bn}), hoping that the new effects
induced by the new resonance terms will not change considerably the original and phenomenologically successful ``fit 10''.
But generally this is not guaranteed.

The detailed analysis based on the current experimental data which would take into account simultaneously and consistently
various resonance estimates for both $O(p^4)$ and $O(p^6)$ LECs is beyond the scope of this paper (however such a project is under investigation \cite{Bijnens:2009hy,Bijnens:2011tb}).
Instead we present another representation that can be used for analyzing the data without addressing the values of the individual $C_i$s.

\section{Dispersive construction}\label{sec:dispersive construction}
The dispersive construction to be presented below is based on the reconstruction theorem
\cite{Stern:1993rg,Knecht:1995tr,Zdrahal:2008bd}, which takes into account only the most general properties of the
amplitude, namely, relativistic invariance, unitarity, analyticity and crossing, supplied with chiral counting
(e.g.\ expansion in powers of momenta and of masses of the pseudoscalars). This framework provides the most general form of the amplitude under consideration in the low-energy domain, up to a remainder of the chiral order $O(p^{8})$. Such a construction requires at the same time the scattering amplitudes related to the original one by two-particle unitarity and by crossing. (Contributions to the unitarity condition arising from intermediate states with more than two pseudo-Goldstone particles only start at $O(p^8)$ --- cf.\ \cite{Stern:1993rg,Knecht:1995tr,Zdrahal:2008bd}.) These amplitudes are constructed along the same lines.
The details of the construction, including a full isospin breaking arising from $m_{\pi^0}\neqm_{\pi^\pm}$, will be published elsewhere \cite{kknzprep} (cf.~also \cite{nas-proceeding2}). In this work we
concentrate on the qualitative description of the result.

The dispersively constructed scattering amplitudes of the pseudo-Goldstone bosons (pGB) take the following general form
\begin{equation}
\mathcal{A}(s,t,u)=\mathcal{N}\left( \mathcal{P}(s,t,u)+\mathcal{U}(s,t,u)\right) +O(p^{8}).  \label{tvar amplitudy}
\end{equation}
Here $\mathcal{N}$ is an overall normalization and $\mathcal{P}(s,t,u)$ is a
third order polynomial with the same symmetry properties with respect to $s$, $t$ and $u$ as the complete amplitude $\mathcal{A}(s,t,u)$. The coefficients of these polynomials for the independent amplitudes related by
two-particle unitarity in all the crossed channels are identical for all the amplitudes and are the only free
parameters entering the game. The non-analytic unitarity part $\mathcal{U}(s,t,u)$, which takes into account the contribution of the two-particle intermediate states in all the crossed channels, is then a known function of these parameters. In the low-energy region, intermediate states containing more than two pGB states contribute only to the $O(p^8)$ remainder, while intermediate states involving other hadronic states contribute to the coefficients of the subtraction polynomial.

In the case of the amplitudes concerning one $\eta$ and three pion states, there are several two-pGB intermediate states to consider: $\pi\pi$, $KK$, $\pi\eta$. Since we shall only be concerned by the decay region, only the nearest singularity, coming from the cut produced by the $\pi\pi$ intermediate state, will be close enough to affect sizably the amplitude. The contribution from the remaining states ($KK$, $\pi\eta$) can be expanded in a polynomial, which is included in ${\cal P}(s,t,u)$ (see also the discussion at the beginning of the next section). Of course, such an approximation would not be appropriate\footnote{However, the presented construction can be extended also to include the unitarity cuts from the other two-pGB intermediate states which are relevant in the scattering region \cite{Zdrahal:2008bd}. This then, however, brings into a game more free parameters (describing such intermediate processes).} to describe the $\eta\pi\to\pi\pi$ amplitude in the scattering region. In conclusion, for our purposes the only relevant related amplitude is therefore the $\pi \pi $ scattering one.

For the charged $\eta \rightarrow 3\pi $ decay channel the polynomial $\mathcal{P}_x(s,t,u)$ can be expressed in terms of six free parameters corresponding to the $t-u$ symmetric expansion at the center of the Dalitz plot
\begin{equation}
\mathcal{P}_x(s,t,u)=A_x m_\eta^{2}+B_x(s-s^{c})+C_x(s-s^{c})^{2}
+D_x \left[(t-s^{c})^{2}+(u-s^{c})^{2}\right]+E_x (s-s^{c})^{3}
+F_x\left[(t-s^{c})^{3}+(u-s^{c})^{3}\right],
\label{polstu}
\end{equation}
which is closely related to the traditional PDG parametrization of the Dalitz plot distribution. We take the overall normalization as
\begin{equation}
{\cal N}_\eta = \frac{\sqrt3}{4}\frac{1}{R} \frac{1}{F_\pi^2}\,,
\end{equation}
so we have simply (cf. (\ref{ARM}))
\begin{equation}\label{vyjadreni M}
{\cal M}_x(s,t,u) = {\cal P}_x(s,t,u) + {\cal U}_x(s,t,u).
\end{equation}
Let us make one remark concerning the Dalitz plot parametrization. Between the polynomial~(\ref{polstu}) and the \emph{linear} parametrization~(\ref{lindalnab}) there is a simple connection. However, the dependence of $\bar a$, $\bar b$, $\bar d$, $\bar f$, $\bar g$ on parameters $A_x$, $B_x$, $C_x$ and $D_x$ is complicated by the presence of these four parameters also in the unitarity part ${\cal U}_x(s,t,u)$ (see below). The direct correspondence can be, however, established for the dependence of $\bar f$, $\bar g$ on $E_x$ and $F_x$ with very simple form
\begin{equation}
\bar f- \bar g \sim E_x + 2 F_x,\label{Eplus2F}
\end{equation}
which we will need in Sec.~\ref{sec:ObOC} (the exact connection will not be needed).

For the related
$\pi ^{+}\pi ^{-}$ $\rightarrow \pi ^{0}\pi ^{0}$ scattering amplitude (which is the only independent one in
the isospin conservation case) we choose the following parametrization of
the polynomial part in terms of the subthreshold parameters \cite{Stern:1993rg,Knecht:1995tr}
\begin{multline}
\mathcal{P}_{\pi \pi }(s,t,u) =\frac{1}{3}\alpha_\pi M_{\pi }^{2}+\beta_\pi \left(s-\frac{4}{3}M_{\pi }^{2}\right) +\frac{\lambda_{1}}{F_{\pi }^{2}}\left(s-2M_{\pi }^{2}\right)^{2}
+\frac{\lambda _{2}}{F_{\pi }^{2}}\left[ \left(t-2M_{\pi
}^{2}\right)^{2}+\left(u-2M_{\pi }^{2}\right)^{2}\right]\\
+\frac{\lambda _{3}}{F_{\pi }^{4}}\left(s-2M_{\pi }^{2}\right)^{3}
+\frac{\lambda _{4}}{F_{\pi }^{4}}\left[\left(t-2M_{\pi }^{2}\right)^{3}+\left(u-2M_{\pi }^{2}\right)^{3}\right]
\end{multline}
and the overall normalization $\mathcal{N}_{\pi \pi }=F_{\pi}^{-2}$.
The unitarity part of the $\eta \rightarrow 3\pi $ decay amplitude $\mathcal{U}_x$ is then a function of a subset of the above polynomial parameters, namely
\begin{equation}
\mathcal{U}_x=\mathcal{U}_x\left(A_x,B_x,C_x,D_x;\alpha_\pi ,\beta_\pi ,\lambda _{1},\lambda _{2}\right).
\end{equation}
The general form of $\mathcal{U}(s,t,u)$ for the process $\text{AB}\rightarrow \text{CD}$ reads
\begin{equation}
\mathcal{U}_{\text{AB}\rightarrow \text{CD}}(s,t,u)
=W_{S}^{0}(s)+W_{T}^{0}(t)+W_{U}^{0}(u)
+(t-u)W_{S}^{1}(s)+(s-u)W_{T}^{1}(t)+(t-s)W_{U}^{1}(u),
\end{equation}
where the discontinuities of the functions $W_{S,T,U}^{0,1}(s)$ are given in terms of the right-hand cut discontinuities of the $S$ and the $P$ partial waves
$S_{\ell}$, $T_{\ell}$ and $U_{\ell}$ $(\ell=0,1)$ of the processes in the $s$-, the $t$- and the $u$-channels, respectively, as\footnote{Note that in our case we need to continue these discontinuities analytically and they become complex (cf.\ \cite{kknzprep,BronzanKacser,Kacser,KWW}).}
\begin{align}
\disc W_{S}^{0}(s) &= 16\pi\, \Bigg(\disc S_{0}(s)
+\Delta_{\text{AB}}\Delta _{\text{CD}} \disc \frac{S_{1}(s)}{\lambda _{\text{AB}}^{1/2}(s)\lambda
_{\text{CD}}^{1/2}(s)}\Bigg),
\\
\disc W_{S}^{1}(s)&=48\pi s \disc \frac{S_{1}(x)}{\lambda
_{\text{AB}}^{1/2}(x)\lambda _{\text{CD}}^{1/2}(x)}\,,
\end{align}
where $\Delta_{ij}$ and $\lambda_{ij}(s)$ were defined in (\ref{lambda}). Similar relations for $W_{T,U}^{0,1}(s)$ can be obtained by an appropriate permutations of $\text{A},\ldots ,\text{D}$. The right hand cut discontinuities  are iteratively constructed from the generalized two-particle partial-wave
unitarity relations  as described in \cite{Stern:1993rg,Knecht:1995tr,Zdrahal:2008bd}.  The functions
$W_{S,T,U}^{0,1}(s)$ are then reconstructed by means of
appropriately subtracted dispersion relation. Note that such a subtraction prescription is an indivisible part of the definition
of the polynomial part of the amplitude. The first iteration
reconstruct the  amplitude at $O(p^4)$ while the second one yields
the $O(p^6)$ results.

For the $\eta \rightarrow \pi ^{+}\pi ^{-}\pi ^{0}$ decay the above general
form simplifies since there are only two independent masses in the problem
and the amplitude is $t-u$ symmetric. We get
\begin{equation}
\mathcal{U}_x(s,t,u) =W_{S}^{0}(s)+W_{T}^{0}(t)+W_{T}^{0}(u)
 +(s-u)W_{T}^{1}(t)+(s-t)W_{T}^{1}(u),
\label{UetaW}
\end{equation}
where the subscripts $S$, $T$ refer to the $\eta \pi ^{0}\rightarrow
\pi ^{+}\pi ^{-}$ and the $\eta \pi ^{+}$ $\rightarrow \pi ^{+}\pi
^{0}$ channels, respectively.
The relevant discontinuities can be rewritten schematically as
\begin{equation}
\disc W_S^0(s) =\sum_{i=1}^{5}\mathcal{F}_{i}(s)\sum_{k=-1}^{4}a_{k}^{(i)}s^{k}
\end{equation}
and similarly for $W_T^0$ (with coefficients $\bar a_k^{(i)}$), while
\begin{equation}
\disc W_T^1(s)
=\sum_{i=1}^{5}\mathcal{F}_{i}(s)\Bigg(
\sum_{k=-1}^{2}\bar b_{k}^{(i)}s^{k}
+\frac{1}{\lambda (s)}\sum_{k=-1}^{5}\bar c_{k}^{(i)}s^{k}+\frac{1}{x\sigma (s)^{2}}
\sum_{k=0}^{3}\bar d_{k}^{(i)}s^{k}\Bigg).
\end{equation}
Here $a_{k}^{(i)}, \bar a_{k}^{(i)}, \ldots , \bar d_{k}^{(i)}$ are known polynomials of the parameters $\{A_x,B_x,C_x,D_x;\alpha_\pi ,\beta_\pi ,\lambda_{1},\lambda _{2}\}$ and the masses $m_\pi$, $m_\eta$;
$\lambda (s)=\lambda_{\pi \eta }(s)$, $\sigma (s)=\lambda _{\pi \pi }(s)/s$ and $\mathcal{F}_{i}(s)$ represents a set of elementary functions listed in Appendix A.

The corresponding functions $W_{S,T,U}^{0,1}(s)$ are now expressed
in terms of the dispersion integrals (the Hilbert transforms)
$\mathcal{G}_{i}(s)$ of these functions, i.e.
\begin{equation}
\mathcal{G}_{i}(s)=\frac{s^{k_{i}}}{\pi }\int_{4m_{\pi }^{2}}^{\infty }\frac{dx}{x^{k_{i}}}\frac{\mathcal{F}_{i}(x)}{x-s}
\end{equation}
with a suitable number $k_{i}$ of subtractions, $0\leq k_{i}\leq 1$.
The $S$-wave contributions in the $s$- and the $t$-channels are given by
\begin{align}
W_{S}^{0}(s)&=\sum_{i=1}^{5}\mathcal{G}_{i}(s) \sum_{k=-1}^{4}a_{k}^{(i)}s^{k},\\
\qquad W_{T}^{0}(s)&=\sum_{i=1}^{5}\mathcal{G}_{i}(s) \sum_{k=-1}^{4}\bar{a}_{k}^{(i)}s^{k}.
\label{WS0}
\end{align}
The $P$-wave contribution in the $t$-channel is more complicated,
\begin{multline}
W_{T}^{1}(s) =\sum_{i=1}^{5}\biggl\{ \mathcal{G}_{i}(s)
\sum_{k=-1}^{2}b_{k}^{(i)}s^{k}+\mathcal{G}_{i}^{(\lambda
)}(s)\sum_{k=-1}^{3}c_{k}^{(i)}s^{k}
+\mathcal{G}_{i}^{(\sigma
)}(s)\sum_{k=0}^{3}d_{k}^{(i)}s^{k}\biggr\}\\
+\sum_{i=1}^{5}\overline{\mathcal{G}}_{i}^{(\lambda
)}(s)~c_{-1}^{(i)}s^{-1}
+\sum_{i=1}^{5}\widetilde{\mathcal{G}}_{i}^{(\lambda
)}(s)~c_{4}^{(i)}s^{4}+\sum_{i=1}^{5}\widehat{\mathcal{G}}_{i}^{(\lambda
)}(s)~c_{5}^{(i)}s^{5},\label{WT1}
\end{multline}
where
(in the following formulae $m_{\pm }=m_\eta\pm m_\pi$)
\begin{align}\label{derGi-first}
\mathcal{G}_{i}^{(\lambda)}(s) &= \frac{1}{m_{+}^{2}-m_{-}^{2}}
\bigg( \frac{\mathcal{G}_{i}(s)-\mathcal{G}_{i}(m_{+}^{2})}{s-m_{+}^{2}}
-\frac{\mathcal{G}_{i}(s)-\mathcal{G}_{i}(m_{-}^{2})}{s-m_{-}^{2}}\bigg),
\\
\mathcal{G}_{i}^{(\sigma )}(s) &= \frac{\mathcal{G}_{i}(s)-\mathcal{G}_{i}(4M_{\pi }^{2})}{s-4M_{\pi }^{2}},\\
\overline{\mathcal{G}}_{i}^{(\lambda )}(s) &=\mathcal{G}_{i}^{(\lambda)}(s)-\mathcal{G}_{i}^{(\lambda )}(0), \\
\widetilde{\mathcal{G}}_{i}^{(\lambda )}(s) &=\frac{1}{s}\frac{1}{m_{+}^{2}-m_{-}^{2}}
\bigg(m_{+}^{2}\frac{\mathcal{G}_{i}(s)-\mathcal{G}_{i}(m_{+}^{2})}{s-m_{+}^{2}}
-m_{-}^{2}
\frac{\mathcal{G}_{i}(s)-\mathcal{G}_{i}(m_{-}^{2})}{s-m_{-}^{2}}\bigg),
\\
\widehat{\mathcal{G}}_{i}^{(\lambda )}(s) &=\frac{1}{s^{2}}\Biggl( \mathcal{G}_{i}(s)+\frac{1}{m_{+}^{2}-m_{-}^{2}}
\biggl( m_{+}^{4}\frac{\mathcal{G}_{i}(s)-\mathcal{G}_{i}(m_{+}^{2})}{s-m_{+}^{2}}-m_{-}^{4}
\frac{\mathcal{G}_{i}(s)-\mathcal{G}_{i}(m_{-}^{2})}{s-m_{-}^{2}}\biggr) \Biggr).
\label{derGi-last}
\end{align}
The dependence of these functions on $\mathcal{G}_i(s)$ ensures the correct discontinuity of the function $W_T^1(s)$ and in addition is dictated by the requirement that the appropriate behavior in the chiral limit \cite{Knecht:1995tr,kknzprep} is reproduced.

The explicit form of the functions $\mathcal{G}_{i}(s)$ as well as the properties of the Hilbert transform are discussed in Appendices A and B. Here we only illustrate the above general procedure by means of the explicit result of the first iteration corresponding to the $O(p^4)$ part of the amplitude, and briefly discuss the $O(p^6)$ result.

\subsection{\texorpdfstring{$\eta\to3\pi$}{eta -> 3 pi} at one-loop order}
\label{sec:etaoneloop}
At the one-loop order our dispersive representation (\ref{vyjadreni M}) of the amplitude ${\cal M}_x(s,t,u)$ simplifies substantially. The polynomial $\mathcal{P}_x(s,t,u)$ is only of the second order,
\begin{equation}
\mathcal{P}_x(s,t,u)=A_x m_\eta^{2}+B_x(s-s^{c})+C_x(s-s^{c})^{2}
+D_x \left[(t-s^{c})^{2}+(u-s^{c})^{2}\right],
\end{equation}
and of all the functions ${\cal F}_i(s)$ and their Hilbert transforms ${\cal G}_i(s)$ that were introduced in the previous section only the case $i=1$ occurs in the unitarity part (\ref{UetaW}). Besides we only need the first term from (\ref{WT1}).

The single function appearing at $O(p^4)$ is thus
\begin{equation}
\mathcal{G}_1(s)=\frac{s}{\pi}\int_{4M_{\pi }^{2}}^{\infty }\frac{dx}{x}\frac{\sigma(x)}{x-s}=\frac{1}{\pi}\left(2+\sigma(s)\log\frac{\sigma(s)-1}{\sigma(s)+1}\right).
\end{equation}
The form of ${\cal F}_1$ was chosen in order to ensure the relation ${\cal G}_1 = (16 \pi)\bar J(s)$ \cite{Knecht:1995tr} (known also as Chew-Mandelstam function \cite{Chew:1960iv}).

The form of the unitarity part (\ref{UetaW}) at the $O(p^4)$ order is extremely simple in this formalism. For the polynomials introduced in~(\ref{WS0}) and~(\ref{WT1}) in the case of the charged decay $\eta\to \pi^+\pi^-\pi^0$ we find
\begin{align}
16\pi F_\pi^2 a^{(1)}_0  &= \frac16 A_x (7 \alpha_\pi - 16 \beta_\pi) m_\pi^2 m_\eta^2
- \frac2 9B_x (\alpha_\pi-\beta_\pi) m_\pi^2 ( 3m_\pi^2 + m_\eta^2),
\label{a01op4}\\
16\pi F_\pi^2 a^{(1)}_1  &= 2 A_x \beta_\pi m_\eta^2 + \frac16 B_x \bigl(4 \alpha_\pi m_\pi^2
- \beta_\pi (7 m_\pi^2 +m_\eta^2)\bigr),
\\
16\pi F_\pi^2 a^{(1)}_2 &= \frac12 B_x \beta_\pi
\end{align}
for the polynomials of S-wave in $s$ channel (the polynomials that are not displayed here are identically zero). And then similarly for S-wave in $t$-channel
\begin{align}
16\pi F_\pi^2\bar a^{(1)}_0  &= \frac13 A_x (\alpha_\pi + 2 \beta_\pi) m_\pi^2 m_\eta^2
 + \frac{1}{18}B_x (\alpha_\pi+2\beta_\pi) m_\pi^2 ( 3m_\pi^2 + m_\eta^2),
\\
16\pi F_\pi^2\bar a^{(1)}_1  &= -\frac12 A_x \beta_\pi m_\eta^2
- \frac{1}{12} B_x \bigl(2 \alpha_\pi m_\pi^2 + \beta_\pi (7 m_\pi^2 +m_\eta^2)\bigr),
\\
16\pi F_\pi^2\bar a^{(1)}_2 &= \frac14 B_x \beta_\pi
\end{align}
and finally the polynomials of the P-wave contributions that are not zero are given by
\begin{align}
16\pi F_\pi^2 \bar{b}^{(1)}_0  &= - \frac13 B_x \beta_\pi m_\pi^2,\\
16\pi F_\pi^2 \bar{b}^{(1)}_1  &= - \frac{1}{12} B_x \beta_\pi.
\end{align}

\subsection{\texorpdfstring{$\eta\to3\pi$}{eta -> 3 pi} at two-loop order}
The whole amplitude at two loops, or equivalently at $O(p^6)$ order, is of course more complicated. We employ here the full form of the polynomial (\ref{polstu}). The non-trivial part follows from the same general form (\ref{UetaW}) with the functions $W_{S,T,U}^{0,1}(s)$ from (\ref{WS0})--(\ref{WT1}), but contrary to the one-loop situation, where we have only one function ${\cal G}_1(s)$, we have to deal with five basic functions ${\cal G}_i$, together with five derived types~(\ref{derGi-first})--(\ref{derGi-last}).
Let us explicitly write down the first coefficient (which stands in front of ${\cal G}_1$ in the $s$-channel of S partial wave and thus together with~(\ref{a01op4}) represents the full $a_0^{(1)}$ at $O(p^6)$):
\begin{widetext}\begin{equation}
\begin{split}
(16\pi)^3&\Delta a_0^{(1)} = \frac{4m_\eta^2}{27F_\pi^4}A_xm_\pi^2
\bigl[
 \left(1152 \alpha_\pi ^2-2619 \alpha_\pi  \beta_\pi +3130 \beta_\pi ^2\right)m_\pi^2-45 \beta_\pi  (\alpha_\pi
   +2 \beta_\pi )m_\eta^2
   +1152 \pi ^2
   \left(17 \lambda _1+18 \lambda _2\right) m_\pi^2\bigr]\\
&+\frac{4}{81F_\pi^4} B_x m_\pi^2 \bigl[
\beta_\pi
    (61 \alpha_\pi -514 \beta_\pi )m_\eta^4+ \left(-459 \alpha_\pi
   ^2+1170 \alpha_\pi  \beta_\pi -1996 \beta_\pi ^2\right)m_\pi^2 m_\eta^2\\
&\qquad
-3 \left(447 \alpha_\pi ^2-987
   \alpha_\pi  \beta_\pi +2150 \beta_\pi ^2\right) m_\pi^4
   -4608 \pi ^2 \left(2 \lambda _1+3
   \lambda _2\right) m_\pi^2 \left(3 m_\pi^2+m_\eta^2\right)\bigr]\\
&+\frac{128}{27F_\pi^2} \pi ^2 C_x
   m_\pi^2 \bigl[7 \alpha_\pi  \left(3 m_\pi^2+m_\eta^2\right)^2-2 \beta_\pi
   \left(81 m_\pi^4+30 m_\pi^2 m_\eta^2+17
   m_\eta^4\right)\bigr]\\
&+\frac{256}{27F_\pi^2} \pi ^2 D_xm_\pi^2 \bigl[7
   \alpha_\pi  \left(3 m_\pi^2+m_\eta^2\right)^2-\beta_\pi  \left(171 m_\pi^4+42
   m_\pi^2 m_\eta^2+43 m_\eta^4\right)\bigr].
\end{split}
\end{equation}\end{widetext}
From this example one can infer the general structure of all other parameters $a_{k}^{(i)}, \ldots , \bar d_{k}^{(i)}$. The full form can be obtained from the authors upon request.

\section{Connection with ChPT: order-by-order correspondence}
\label{sec:conwithchpt}
Let us briefly comment on the connection of the dispersive construction with
the standard ChPT expansion. In analogy to the dispersive one, the $O(p^{6})$ ChPT
amplitude can also be split into a polynomial part and a non-analytic
unitarity part. The former corresponds to the tree-level counterterm
contributions as well as to the chiral logs and sunset graphs, while the
latter takes explicitly into account the nontrivial contributions of the
loops. Though this splitting is not unambiguous and depends on the particular
definition of the nontrivial part of the loop graphs, the unitarity part
has to reproduce the correct discontinuities of the amplitude as required by
(generalized) unitarity and corresponding to the two-particle intermediate
states. Along with the pure pion loop contributions also the higher
intermediate states are taken into account, namely, the graphs with kaons and
$\eta $ inside the loops. However, below the $\pi \eta $ threshold the
contributions of discontinuities corresponding to the $\pi \eta $, $KK$ and $\eta \eta $ intermediate states are analytic and can therefore be expanded in powers of $s$, $t$, $u$. Sufficiently far below these thresholds one can show that their effects can be approximated by means of only the terms up to the third order (cf.\ \cite{nas-proceeding3} and the numerical estimate of such error made in Sec.~\ref{sec:ObOC}). As a result we should obtain in this region an
approximate ChPT amplitude with the same structure as our dispersively
constructed amplitude (recall that both of them include the higher non-Goldstone
intermediate states contributions only effectively through the low-energy and the subtraction constants, respectively). The only difference is that the polynomial part of the $O(p^{6})$ ChPT amplitude is generally
complex due to the contribution of the sunset diagram\footnote{Note that this diagram ($\eta\to3\pi\to3\pi$) does not contribute to the unitarity cut of the $\eta\pi\to\pi\pi$ amplitude but instead its contribution in the decay region is analytic and can be expanded into polynomial. This polynomial can be complex since $m_\eta$ is unstable ($m_\eta>3m_\pi$).}  with three
intermediate pions which develop nonzero imaginary part. However, it has been
found to be tiny in \cite{BG,nas-proceeding3} and therefore can be neglected. We reverify this observation in Sec.~\ref{sec:ObOC}.

These common features of both amplitudes suggest that the $O(p^{6})$ ChPT
amplitude ${\mathcal{A}}_{\mathrm{ChPT}}(s,t,u)$, which we write in the form
\begin{equation}
{{\mathcal{A}}}_{\mathrm{ChPT}}(s,t,u)={\cal N}_\eta\bigg( \mathcal{M}_{\mathrm{ChPT}}^{(2)}(s,t,u)
+\mathcal{M}_{\mathrm{ChPT}}^{(4)}(s,t,u)+\mathcal{M}_{\mathrm{ChPT}}^{(6)}(s,t,u)\bigg) ,
\end{equation}
can be reproduced as a special case of the dispersively constructed one. This can be quantified as follows in terms of what we call \emph{order-by-order fit}. The ChPT amplitude in our dispersive parametrization is then represented by expressing particular chiral orders of our subtraction constants $A_x, \dots, F_x$ and $\alpha_\pi, \dots, \lambda_{2}$ in terms of the LECs of ChPT, quark masses and chiral logarithms. Such expressions are then useful when one wants to organize the chiral result and to identify the renormalization-scale invariant combinations of LECs on which the amplitude depends. For the aims of the current work, it is however sufficient to perform this matching numerically and obtain the numerical values of our subtraction constants using the procedure described in the following lines (note that the same procedure would remain valid also if we wanted to obtain the analytic expressions, but instead of fitting the numerical results we would just compare expressions coming from ChPT with the ones of the analytic dispersive construction).

Let us formally split the parameters $A_x$, \ldots ,$F_x$
of our amplitude into their $O(p^{2})$, $O(p^{4})$ and $O(p^{6})$ parts, i.e.
\begin{align}
A_x &=A_x^{(2)}+\Delta A_x^{(4)}+\Delta A_x^{(6)},\\
B_x &=B_x^{(2)}+\Delta B_x^{(4)}+\Delta B_x^{(6)},\\
C_x &=C_x^{(4)}+\Delta C_x^{(6)},\\
D_x &=D_x^{(4)}+\Delta D_x^{(6)},\\
E_x &\equiv E_x^{(6)},\\
F_x &\equiv F_x^{(6)}.
\end{align}
This induces a following splitting of the polynomial part of the
amplitude
\begin{equation}
{\mathcal{P}}_x(s,t,u)={\mathcal{P}}_x^{(2)}(s,t,u)+{\mathcal{P}}_x^{(4)}(s,t,u)+{\mathcal{P}}_x^{(6)}(s,t,u),
\end{equation}
where
\begin{align}
{\mathcal{P}}_x^{(2)}(s,t,u) &=A_x^{(2)}m_\eta^{2}+B_x^{(2)}(s-s^{c}), \\
{\mathcal{P}}_x^{(4)}(s,t,u) &=\Delta A_x^{(4)}m_\eta^{2}+\Delta
B_x^{(4)}(s-s^{c})+C_x^{(4)}(s-s^{c})^{2}+D_x^{(4)}\left[(t-s^{c})^{2}+(u-s^{c})^{2}\right],  \\
\begin{split}
{\mathcal{P}}_x^{(6)}(s,t,u) &=\Delta A_x^{(6)}m_\eta^{2}+\Delta
B_x^{(6)}(s-s^{c})
+\Delta C_x^{(6)}(s-s^{c})^{2}
+\Delta D_x^{(6)}\left[ (t-s^{c})^{2}+(u-s^{c})^{2}\right]\\
&+E_x^{(6)}(s-s^{c})^{3}+F_x^{(6)}\left[(t-s^{c})^{3}+(u-s^{c})^{3}\right].
\end{split}\label{P6}
\end{align}
Note that the unitarity part ${\mathcal{U}}_x(s,t,u)$ splits by construction naturally into the genuine one-loop $O(p^{4})$ and the remaining $O(p^{6})$ parts that correspond to the first and the second iteration of the generalized
unitarity relations, respectively, (see \cite{kknzprep} for more details),
\begin{equation}
{\mathcal{U}}_x={\mathcal{U}}_x^{(4)}\left(A_x,B_x;\alpha_\pi ,\beta_\pi \right)
+\mathcal{U}_x^{(6)}\left(A_x,B_x,C_x,D_x;\alpha_\pi ,\beta_\pi
,\lambda _{1},\lambda _{2}\right).
\end{equation}
The unitarity part ${\mathcal{U}}_x^{(4)}$ has been given in Sec.~\ref{sec:etaoneloop}, where we have written out the explicit dependence on the polynomial parameters of the $\eta \rightarrow 3\pi $ and $\pi \pi
\rightarrow \pi \pi $ amplitudes. The $O(p^{6})$ part consists further of the
genuine two-loop part and the one-loop part
\begin{equation}
{\mathcal{U}}_x^{(6)}(A_x,B_x,C_x,D_x;\alpha_\pi ,\beta_\pi
,\lambda _{1},\lambda _{2})
={\mathcal{U}}_{\mathrm{2-loop}}^{(6)}(A_x,B_x;\alpha_\pi ,\beta_\pi )
+{\mathcal{U}}_{\mathrm{1-loop}}^{(6)}(A_x,B_x,C_x,D_x;\alpha_\pi ,\beta_\pi ,\lambda _{1},\lambda _{2}).
\end{equation}

The $O(p^{2})$ ChPT amplitude ${\mathcal{M}}_{\mathrm{ChPT}}^{(2)}(s,t,u)$
is now exactly reproduced by ${\mathcal{P}}_x^{(2)}(s,t,u)$ with
\begin{equation}
A_x^{(2)} =\frac{m_\eta^{2}-m_{\pi }^{2}}{3m_\eta^{2}},\qquad B_x^{(2)} =1.\label{ABLO}
\end{equation}
The imaginary part of the $O(p^{4})$ ChPT amplitude below the $\pi \eta $
threshold is fixed by unitarity and therefore there holds exactly
\begin{equation}
\im {\mathcal{M}}_{\mathrm{ChPT}}^{(4)}(s,t,u)=\im \mathcal{U
}_x^{(4)}(A_x^{(2)},B_x^{(2)};\alpha_\pi ^{(2)},\beta_\pi ^{(2)}),
\end{equation}
where
\begin{equation}
\alpha_\pi ^{(2)}=\beta_\pi ^{(2)}=1 \label{abLO}
\end{equation}
are the leading order ChPT values of the $\pi \pi \rightarrow \pi \pi $
subthreshold parameters. Hence, up to a polynomial of the second order
in $s$, $t$ and $u$, the amplitudes ${\mathcal{M}}_{\mathrm{ChPT}}^{(4)}(s,t,u)$ and ${\mathcal{U}}_x^{(4)}(A_x^{(2)},B_x^{(2)};\alpha_\pi ^{(2)},\beta_\pi ^{(2)})$ coincide (here we have tacitly assumed
that the higher two-particle intermediate states contributions to $\mathcal{M}_{\mathrm{ChPT}}^{(4)}(s,t,u)$ has been expanded to the second order in $s$, $t$ $\ $and $u$ as described above) and we can therefore write
\begin{equation}
{\mathcal{M}}_{\mathrm{ChPT}}^{(4)}(s,t,u)={\mathcal{P}}_x^{(4)}(s,t,u)+{\mathcal{U}}_x^{(4)}(A_x^{(2)},B_x^{(2)};\alpha_\pi ^{(2)},\beta_\pi
^{(2)})
\end{equation}
for appropriate $\Delta A_x^{(4)}$, $\Delta B_x^{(4)}$, $C_x^{(4)}$ and $D_x^{(4)}$. These parameters are found numerically by fitting the difference
\begin{equation}
\Delta ^{(4)}(s,t,u)={\mathcal{M}}_{\mathrm{ChPT}}^{(4)}(s,t,u)-{\mathcal{U}}_x^{(4)}(A_x^{(2)},B_x^{(2)};\alpha_\pi ^{(2)},\beta_\pi ^{(2)})
\end{equation}
to the second order polynomial ${\mathcal{P}}_x^{(4)}(s,t,u)$. When these $O(p^{4})$ parameters are fixed, we proceed similarly to the $O(p^{6})$
order. We compute the $O(p^6)$ corrections to the unitarity part,
\begin{equation}\begin{split}
{\mathcal{V}}_x^{(6)}(s,t,u) &={\mathcal{U}}_x^{(4)}(\Delta A_x^{(4)},\Delta B_x^{(4)};\alpha_\pi ^{(2)},\beta_\pi ^{(2)})
+{\mathcal{U}}_x^{(4)}(A_x^{(2)},B_x^{(2)};\Delta \alpha_\pi ^{(4)},\Delta \beta_\pi
^{(4)})\\
&+{\mathcal{U}}_{\mathrm{1-loop}}^{(6)}(A_x^{(2)},B_x^{(2)},C_x^{(4)},D_x^{(4)};\alpha_\pi ^{(2)},\beta_\pi ^{(2)},\lambda
_{1}^{(4)},\lambda _{2}^{(4)})
+{\mathcal{U}}_{\mathrm{2-loop}}^{(6)}(A_x^{(2)},B_x^{(2)};\alpha_\pi ^{(2)},\beta_\pi ^{(2)}),
\end{split}\end{equation}
where in addition to the parameters known from the previous steps there appear the NLO corrections of the subthreshold parameters of $\pi\pi\to\pi\pi$ scattering that are needed as inputs to this procedure. The discontinuities originating from the $\pi \pi $ intermediate states in \mbox{$s$-}, $t$- and $u$-channels of ${\mathcal{M}}_{\mathrm{ChPT}}^{(6)}(s,t,u)$ and of this ${\mathcal{V}}^{(6)}_x(s,t,u)$ coincide (modulo a power expansion of the higher-intermediate-state
contributions to the third power as discussed above). Finally, we fit the difference
\begin{equation}
\Delta ^{(6)}(s,t,u)={\mathcal{M}}_{\mathrm{ChPT}}^{(6)}(s,t,u)-{\mathcal{V}}_x^{(6)}(s,t,u)
\end{equation}
to the third order polynomial ${\mathcal{P}}_x^{(6)}(s,t,u)$ and set the remaining $O(p^6)$ parameters. In this way, all
the parameters of the polynomial part of the amplitude are numerically
determined and the $O(p^{6})$ ChPT amplitude ${\mathcal{A}}_{\mathrm{ChPT}}(s,t,u)$ is represented now as ${\mathcal{A}}_{\mathrm{ChPT}}(s,t,u) \rightarrow {\mathcal{A}}_{\mathrm{ord}}^{\mathrm{disp}}(s,t,u)$, where
\begin{equation}
{\mathcal{A}}_{\mathrm{ord}}^{\mathrm{disp}}(s,t,u)={\mathcal{P}}_x^{(2)}(s,t,u)+{\mathcal{P}}_x^{(4)}(s,t,u)+{\mathcal{P}}_x^{(6)}(s,t,u)
+{\mathcal{U}}_x^{(4)}(A_x^{(2)},B_x^{(2)};\alpha_\pi ^{(2)},\beta_\pi
^{(2)})+{\mathcal{V}}^{(6)}(s,t,u).
\end{equation}

By construction, the chiral orders of the various contributions to $\mathcal{A}_{\mathrm{ChPT}}^{\mathrm{disp}}(s,t,u)$ were strictly respected --- for instance the genuine two-loop unitarity corrections depend only on
the leading order parameters $A_x^{(2)},B_x^{(2)}$ and
$\alpha_\pi ^{(2)},\beta_\pi ^{(2)}$. However, the known general form of
the dispersively constructed amplitude ${\mathcal{A}}(s,t,u)$ can
be further used in order to go beyond the strict chiral expansion
and partially resum also the higher chiral-order contributions.
This representation that we call \emph{resummed fit} can be achieved by means of inserting the full parameters
$A_x$, \ldots ,$F_x$ obtained by the above
\emph{order-by-order fit} and the full $O(p^{4})$ $\pi \pi$ subthreshold
parameters (or even the experimental values of the $\pi\pi$ subthreshold parameters from \cite{DescotesGenon:2001tn})
into the unitarity part of the amplitude, i.e.\ to define
\begin{equation}
{\mathcal{A}}_{\mathrm{res}}^{\mathrm{disp}}(s,t,u)={\mathcal{P}}_x^{(2)}(s,t,u)+{\mathcal{P}}_x^{(4)}(s,t,u)+{\mathcal{P}}_x^{(6)}(s,t,u)
+{\mathcal{U}}_x^{(4)}(A_x,B_x;\alpha_\pi ,\beta_\pi )
+{\mathcal{U}}_x^{(6)}(A_x,B_x,C_x,D_x;\alpha_\pi ,\beta_\pi ,\lambda
_{1},\lambda _{2}).
\end{equation}
The difference ${\mathcal{A}}_{\mathrm{res}}^{\mathrm{disp}}-{\mathcal{A}}_{\mathrm{ord}}^{\mathrm{disp}}$ is of order $O(p^{8})$ and contains effectively
contributions of the one and the two-loop graphs with higher-order counterterms.
It might be therefore treated as a rough estimate of the convergence of the
chiral expansion.

Let us note that we could also use another parametrization of the relevant $\pi \pi $ scattering amplitude based on the scattering lengths and effective
ranges instead of the subthreshold parameters (see \cite{kknzprep,nas-proceeding2} for details) and
repeat the above construction along the same lines. In such a case the
amplitude ${\mathcal{A}}_{\mathrm{ord}}^{\mathrm{disp}}(s,t,u)$
has to be numerically the same as before, namely, the parameters $A_x$, \ldots ,$F_x$ should be the same. However, the amplitude will now
depend on the scattering lengths and the effective ranges of the $\pi \pi $
scattering taken up to the order $O(p^{4})$. Provided we then use the
experimental values of these parameters in the \emph{resummed} amplitude $\mathcal{A}_{res}^{\mathrm{disp}}(s,t,u)$, we can interpret the result as a partial
resummation of the two-particle rescattering in the final state. The
numerical effect of such a resummation might be even larger than within the
previous parametrization, because the scattering lengths are known to have
much worse convergent chiral expansion than the subthreshold parameters.

\section{Analysis of the charged decay: \texorpdfstring{$\eta\to\pi^+\pi^-\pi^0$}{eta -> pi(+) pi(-) pi(0)}}\label{sec:charged}
We have prepared everything to employ the dispersive representation for our analysis of the process $\eta\to3\pi$.
It proceeds as follows. We start with the NNLO result of ChPT \cite{BG}. We determine the values of our parameters that reproduce the ChPT result, thereby checking also that the correspondence between these two frameworks holds using the \emph{order-by-order fit} as outlined in the previous section. Our further analysis is motivated by the conclusion of Sec.~\ref{sec:Ci} that the observed mismatch between the ChPT $O(p^6)$ predictions of the Dalitz parameters and their experimental determination by KLOE might be caused by the incorrect determination of the $O(p^6)$ LECs $C_i$ of ChPT. We therefore study the dispersive representation of ChPT with the values of the $C_i$s undetermined and try to find the values of their combinations that reproduces the experimental data. Finally, after that we change completely the strategy and fit directly our dispersive representation to the experimental data. Such a fit gives us the $\eta\to3\pi$ amplitude up to the normalization that is determined from the matching with ChPT in the region where we can believe the ChPT result. In all the cases we are interested in the distribution we obtain and then by comparing the decay widths computed from these distributions (by integration of the square of the amplitude over the physical phase space) with the experimentally measured one, we obtain the value of $R$.

In principle, this could be done for both the charged and the neutral $\eta\to3\pi$ decay. However, as was discussed in Sec.~\ref{Dalitz plot parametrization}, no current experiment determined more than just one Dalitz parameter describing the neutral decay, thus we concentrate mainly on the charged one. Even in the charged sector the experimental situation is poor --- only KLOE \cite{KLOE-etac} published just 4+1 Dalitz parameters (the last one claimed to be compatible with zero) describing the amplitude. From these values of the Dalitz parameters we have constructed a distribution in the physical region (in the similar way as done in \cite{Kupsc:2009zza}) and all our experimental fits are fits to such KLOE-like distributions, in our analysis we therefore depend fully on these KLOE measurements.

\subsection{Order-by-order correspondence: obtaining numerical ChPT distribution}
\label{sec:ObOC}
As was discussed in the previous section we can obtain the approximate chiral $O(p^6)$ amplitude as a special case of our dispersive parametrization with some particular values of our parameters. The correspondence between such amplitude and the result of ChPT has to be almost identical neglecting only small effects descending from expansion of the two-kaon and the $\pi\eta$ contributions and a tiny imaginary part produced by sunset-like diagrams. In principle, working in the $(s,t)$ plane they should agree in the region%
\footnote{Note that even though we are talking about the expansion for small Mandelstam variables (e.g.\ $s$ and $t$), it does not simply mean that the smaller these variables are, the better agreement between these theoretical frameworks we obtain. The amplitude depends on three kinematic variables $s,t,u$ which are connected by relation~(\ref{soucet Mandelstamu decay}). So keeping two of them small, the third one is shifted up by the $m_\eta$ mass.}
for small $s,t$ under the $\pi\eta$ thresholds in all the crossed channels. Although much bigger deviation should be visible only after $KK$ threshold (the contribution of $\pi\eta$ is very small) we stick on this as a strict limit of our method.
Influence of systematic uncertainties is studied using different regions in our matching procedure (see below).
The physical and the matching regions together with the $\pi\eta$ threshold are depicted in Fig.~\ref{fig:domains}.

\begin{figure}
\begin{center}
\includegraphics{tabbijst.eps}
\end{center}
\caption{Domain of applicability of the dispersive method together with the
physical region (ellipse). Full magenta line
represents the $\pi\eta$ threshold in the $u$-channel (the only $\pi\eta$ threshold visible for this range of $s$ and $t$). The blue dot-dashed line represents the axis of the $t-u$ symmetry of the amplitude: e.g.\ the Adler zero in $s,t=0$ (represented by a circle) has its counterpart for $s=0$ and $t=m_\eta^2 + 3 m_\pi^2$.
Dotted lines denote the $4 m_\pi^2$ thresholds in the $s$- and the $t$-channels.
Dashed lines define different regions in the matching procedure (see main
text).}
\label{fig:domains}
\end{figure}

We match the amplitude along the lines of the previous section really order by order.
The correspondence for LO and the imaginary part of NLO can be verified analytically, having $A_x^{(2)}$ and $B_x^{(2)}$
from (\ref{ABLO}) and $\alpha_\pi^{(2)}=\beta_\pi^{(2)}=1$. After that we have proceeded with the matching numerically.
From the NLO real part of the amplitudes we have fitted the parameters $A_x^{(4)}$, $B_x^{(4)}$, $C_x^{(4)}$, $D_x^{(4)}$
(in the notation of Sec.~\ref{sec:conwithchpt}, $A_x^{(4)} = A_x^{(2)}+\Delta A_x^{(4)}$, etc.). After that we have verified the matching of the imaginary NNLO amplitudes and finally from the real part of the NNLO amplitudes fitted the parameters
$A_x^{(6)}, \dots, F_x^{(6)}$ (again $A_x^{(6)} = A_x^{(4)}+\Delta A_x^{(6)}$; this superscript is used for the NNLO values of $A_x, \dots$ just to distinguish these values from the ones of the \emph{overall fit} from Sec.~\ref{sec:res}).

Concerning the $\pi\pi$ part we follow closely the determination of its subthreshold parameters as established in~\cite{Knecht:1995tr}. For the particular values we have used~(\ref{abLO}) for the leading order and set the NLO values to be
\begin{equation}\begin{split}
      \alpha_\pi^{(4)} = 1.044\,,&\qquad
      \beta_\pi^{(4)}  = 1.083\,,\\
      \lambda_1^{(4)}  = -1.43\times 10^{-3}\,,&\qquad
      \lambda_2^{(4)}  = 8.5\times 10^{-3}.
\end{split}\end{equation}
The fits were performed in the following regions\footnote{Since the amplitude is $t\leftrightarrow u$ symmetric, one can fit it only in the region below the $t-u$ line.} (all numbers in GeV$^2$; cf.~Figure~\ref{fig:domains}):
\begin{itemize}
 \item set 0: the physical region;
 \item set 1: the square region around the Adler zero $(s,t)\in (-0.05,0.05)$;
 \item set 2: the triangle region between the lines $s=4 m_\pi^2$, $t=4 m_\pi^2$ and the $\pi\eta$ threshold;
 \item set 3: $s\in (-0.05,0.05)$, $t$ between the $\pi\eta$ threshold and the $t\leftrightarrow u$ axis;
 \item set 4: $s\in (-0.1,0.1)$, $t$ between the $\pi\eta$ threshold and the $t\leftrightarrow u$ axis.
\end{itemize}
Distances between the points in grids are constant in both the $s$ and the $t$ directions,
and the approximate total number of them is the following: set 1: 300, set 2: 900, set 3:
1600 and set 4: 4400. Further, for the physical region (set 0) we have chosen the same points that were used
in~\cite{BG}, i.e.\ 174 points, which is a very similar number to the KLOE's
number of bins (154, cf.\ also discussion in \cite{Kupsc:2009zza}).
The different regions with the different numbers of points were set in order to have systematic and statical errors under control.

For the fits we have used \textsf{MINUIT} package with the weights of the individual points set to $\Delta O(p^6)/2$.
Results for the NLO parameters are summarized in Table~\ref{tab:sub4} and the ones for the NNLO parameters then in
Table~\ref{tab:sub6}. The error bars quoted for the individual parameters are results of \textsf{MINUIT}.

\begin{table*}
\centering
\begin{tabular}{@{\hspace{1em}}c@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}
l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{0.5em}}}
\hline\hline
& \multicolumn{2}{c@{\hspace{2em}}}{set 0} & \multicolumn{2}{c@{\hspace{2em}}}{set 1} & \multicolumn{2}{c@{\hspace{2em}}}{set 2} & \multicolumn{2}{c@{\hspace{2em}}}{set 3} & \multicolumn{2}{c@{\hspace{0.5em}}}{set 4}\\\hline
$A_x^{(4)}$
& $0.464$ & $ 0.017$
& $0.457$ & $ 0.040$
& $0.459$ & $0.011$
& $0.452$ & $0.016$
& $0.462$ & $0.010$
\\
$B_x^{(4)}$
& $1.95$ & $0.15$
& $1.90$ & $0.21$
& $1.91$ & $0.05$
& $1.88$ & $0.09$
& $1.93$ & $0.06$
\\
$C_x^{(4)}$
& $-0.42$ & $ 4.4$
& $-0.68$ & $0.91$
& $-0.62$ & $0.17$
& $-0.76$ & $0.41$
& $-0.51$ & $0.25$
\\
$D_x^{(4)}$
& $1.07$ & $ 3.8$
& $1.04 $ & $ 0.03$
& $1.04 $ & $ 0.01$
& $1.04 $ & $ 0.02$
& $1.04 $ & $ 0.01$
\\
\hline\hline
\end{tabular}
\caption{The values of the $O(p^4)$ dispersive parameters (in appropriate powers of [$GeV$]) corresponding to matching to ChPT NLO amplitude.}
\label{tab:sub4}
\end{table*}
\begin{table*}
\centering
\begin{tabular}{@{\hspace{1em}}c@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{0.5em}}}
\hline\hline
& \multicolumn{2}{c@{\hspace{2em}}}{set 0} & \multicolumn{2}{c@{\hspace{2em}}}{set 1} & \multicolumn{2}{c@{\hspace{2em}}}{set 2} & \multicolumn{2}{c@{\hspace{2em}}}{set 3} & \multicolumn{2}{c@{\hspace{0.5em}}}{set 4}\\\hline
$A_x^{(6)}$
& $0.577 $ & $ 0.013$
& $0.581$ & $ 0.003$
& $0.581$ & $ 0.023$
& $0.577$ & $ 0.002$
& $0.583$ & $ 0.011$
\\
$B_x^{(6)}$
& $2.42 $ & $ 0.22$
& $2.460 $ & $ 0.012$
& $2.47$ & $ 0.19 $
& $2.44$ & $ 0.01$
& $2.49$ & $  0.10 $
\\
$C_x^{(6)}$
& $0.24 $ & $ 3.4$
& $0.30 $ & $  0.11$
& $0.38$ & $ 1.7 $
& $0.20 $ & $ 0.09$
& $0.55$ & $  0.90 $
\\
$D_x^{(6)}$
& $1.55 $ & $ 2.4$
& $ 1.57 $ & $ 0.02$
& $1.58$ & $ 0.03$
& $1.58 $ & $ 0.02$
& $1.58$ & $  0.02$
\\
$E_x^{(6)}$
& $5 $ & $  149$
& $5.4 $ & $  0.7$
& $5.6$ & $ 4.7 $
& $5.1 $ & $ 0.6 $
& $6.1$ & $  2.6$
\\
$F_x^{(6)}$
& $-4 $ & $  84$
& $-3.6 $ & $  0.1$
& $-3.7$ & $ 0.2$
& $-3.7 $ & $ 0.1$
& $-3.7$ & $  0.1$
\\\hline\hline
\end{tabular}
\caption{The values of the $O(p^6)$ dispersive parameters (in appropriate powers of [$GeV$]) corresponding to matching to ChPT NNLO amplitude.}
\label{tab:sub6}
\end{table*}

At the moment we have in hands the dispersively constructed amplitude (i.e.\ the analytic formula) which is numerically equivalent (or very close) to NNLO ChPT amplitude. We can verify the equivalence also by computing the decay width.

Our dispersive representation was constructed in accordance with chiral
perturbation theory and we have chosen similar normalization as used in \cite{BG} (with an extra factor $1/F_\pi^2$).
We can thus compare directly a neat amplitude ${\cal M}_x(s,t,u)$ with the isospin-breaking
factor pulled out as defined in (\ref{ARM}).
The result of the integration of the amplitude square over the physical phase space is (cf. (6.7) in \cite{BG}):
\begin{equation}
 \Gamma(\eta\to\pi^+\pi^-\pi^0) = \sin^2 \epsilon \times 2.68 \text{ MeV}, \quad
\text{[ChPT]}
\end{equation}
where we have introduced
\begin{equation}
\sin\epsilon = \frac{\sqrt{3}}{4R}\,.
\end{equation}
Comparing this result with the experimental measurement for the decay rate \cite{PDG} we
arrive at the value which exactly reproduces the one of \cite{BG} (mind the typo in \cite{BG})
\begin{equation}
 R  = 41.3. \quad
\text{[ChPT]}
\end{equation}

We can use result \cite{BG} also for a numerical estimate of the error induced by a few approximations in our parametrization we have made with respect to the ChPT computation. As was discussed in Sec.~\ref{sec:conwithchpt}, we have neglected the imaginary parts of our parameters (which are connected with the contribution of the sunset diagram). In the physical region we have performed fits, in which we have allowed the $O(p^6)$ parameters to be complex. We have found that the NNLO ChPT result is very well approximated by adding a constant imaginary term $\im A_x^{(6)}=0.080\pm0.064$. By neglecting this term in the computation of $R$ we introduce an error of 0.1\%. Similarly, we have neglected higher than third order polynomial terms in the expansion of $KK$ and $\pi\eta$ contributions (in the decay region). We can estimate the corresponding error by addition of some higher-order terms into the polynomial. The symmetries dictate that the fourth order polynomial would contain terms $G_x(s-s^c)^4+H_x(s - s^c)^2 (t - s^c) (u - s^c) + I_x\left((t - s^c)^4 + (u - s^c)^4\right)$. From the dimensional considerations, the contribution of $KK$ intermediate states into these parameters should be $\sim\frac{1}{M_K^4(4\pi F_\pi)^2}\lesssim 10^2\,\mathrm{GeV}^{-6}$ (and similarly for $\pi\eta$), whereas even if all of them were $\sim10^4\,\mathrm{GeV}^{-6}$ the shift in the determined $R$ would be 0.1\%. Both of the errors are therefore negligible with respect to the other sources or error discussed in the following analyses.

\subsection{Correction to order-by-order fit: Correcting the \texorpdfstring{$C_i$s}{Ci's} in ChPT}\label{sec:corOO}
In the previous subsection we have constructed the dispersive amplitude reproducing ChPT in the region where our method is
applicable. It is no surprise that if we fitted this dispersive representation to the Dalitz parametrization (\ref{dalitz etac}) as was done in \cite{BG}, we would obtain the same values of the Dalitz parameters as \cite{BG}. In Sec.~\ref{sec:Ci} we have found an indication that the discrepancy between so obtained values and the values measured by KLOE can be (at least partially) caused by the incorrect values used for the $O(p^6)$ LECs $C_i$ of ChPT. The contribution of the $C_i$s to the $O(p^6)$ amplitude is polynomial and real and so changing them means changing the $O(p^6)$ part of our polynomial (\ref{P6}) --- shifting the parameters appearing in it. By studying the chiral amplitude obtained from our previous analysis with an unknown $O(p^6)$ polynomial added,
\begin{equation}\label{opravena amplitude}
 {\cal A}_x^{cor}(s,t,u) = {\cal A}_x(s,t,u) + \Delta {\cal A}_x(s,t,u),
\end{equation}
with $\Delta {\cal A}_x(s,t,u) = {\cal N}_\eta \Delta {\cal P}_x(s,t,u)$ and
\begin{multline}
 \Delta {\cal P}_x (s,t,u) = \Delta A_x M_\eta^2 + \Delta B_x (s-s^c) +
\Delta C_x(s-s^c)^2
+ \Delta D_x \bigl[(t-s^c)^2 + (u-s^c)^2 \bigr]\\
+ \Delta E_x (s-s^c)^3
 + \Delta F_x \bigl[(t-s^c)^3 + (u-s^c)^3 \bigr]\,,
\label{DPx}
\end{multline}
we can thus study the impact of the corrected $C_i$s on the chiral $O(p^6)$ amplitude.

Provided the dominant part of the discrepancy between the NNLO chiral result and the measured amplitude is hidden just in the incorrect determination of the $C_i$s, the chiral $O(p^6)$ amplitude with the correct set of the $C_i$s, and thereby also the corrected amplitude ${\cal A}_x^{cor}(s,t,u)$, should reproduce the physical data. Therefore by fitting the KLOE-like distribution, we should obtain the values of the dispersive parameters corresponding to the correct values of the $C_i$s. By comparison of these values with the analytic expressions of these parameters in terms of the $C_i$s, one could obtain approximate constraints that the correct values of the $C_i$s should fulfill\footnote{At the current level these constraints could be formulated in terms of reproducing the measured Dalitz plot parameters. For every such parameter by using relations of Section III.A and the observed difference between its experimental value and the value coming from \cite{BG} with all $C_i=0$, one obtains one constraint on the $C_i$s. Note that provided the information on $R$ was supplied from another source with enough accuracy, one could obtain one additional constraint on the $C_i$s. Unfortunately, such constraints are very complicated and would need to be analyzed together with additional constraints coming from other processes (similarly as was done in \cite{Bijnens:2011tb}) in order to provide any useful information on the values of $C_i$s.}.

However, in the case the change of the $C_i$s is insufficient in order to obtain the physically measured amplitude and there still remains a big difference between the physical amplitude and the one obtained from NNLO ChPT, the fit of the amplitude ${\cal A}_x^{cor}(s,t,u)$ to the physical data would mean this time that this difference was parametrized (and approximated) by the polynomial $\Delta {\cal A}_x(s,t,u)$.

Unfortunately, just from the fit of the amplitude ${\cal A}_x^{cor}(s,t,u)$ to the physical data there is no way how to distinguish between these two scenarios --- either that the $O(p^6)$ chiral amplitude with the correct values of the $C_i$ describes well the physical amplitude or that the higher-order remainder can be on the physical region approximated by the polynomial (or if both of the situations are present in a combination, there is no way how to separate these two contributions). We should be therefore careful with the interpretation of the result of such a fit and take this fit just as a starting point for the deeper analysis of the chiral amplitudes. Note that in Sec.~\ref{sec:Ci} we have listed a few criteria that would indicate the realization of the first scenario.

In any case the amplitude constructed that way should describe the data better than the ChPT parametrization of the previous subsection and the determination of $R$ from this distribution will be closer to the real one.

We have therefore fitted all $\Delta$ parameters from (\ref{DPx}) using the condition that the amplitude (\ref{opravena amplitude}) has to fulfill the distribution of data based on KLOE. This distribution is
limited only within the ellipse of the physical region (cf.\ Figure~\ref{fig:domains}).
We have performed the fit for the following two data sets. First we have created exactly the same points as done for set 0 (or equivalently in \cite{BG}) and then much more (2500) in order to study how this affects the dependence on the statistics. Using the different sets of parameters for the chiral amplitude ${\cal A}_x(s,t,u)$ as summarized in Table~\ref{tab:sub4} and \ref{tab:sub6} has very little effect on the resulting $\Delta$ parameters, so we display only their average, cf.\ Table~\ref{tab:cor}. For further applications we have also fitted the complete $O(p^6)$ polynomial
${\cal P}_x^{(6)}(s,t,u)+\Delta {\cal P}_x (s,t,u)$, whose results are presented in the first column of Table~\ref{tab:resummed}.
\begin{table}
\centering
\begin{tabular}{@{\hspace{1em}}c@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{1em}}}\hline\hline
& \multicolumn{2}{c@{\hspace{2em}}}{\# 174} & \multicolumn{2}{c@{\hspace{1em}}}{\# 2500} \\\hline
$\Delta A_x$
& $-0.05 $ & $ 0.3$
& $-0.029$ & $ 0.003$
\\
$\Delta B_x$
& $-0.5 $ & $ 1$
& $-0.46 $ & $ 0.01$
\\
$\Delta C_x$
& $-7 $ & $ 2$
& $-6.97 $ & $ 0.07$
\\
$\Delta D_x$
& $-0.7 $ & $ 0.8$
& $-0.64 $ & $ 0.02$
\\
$\Delta E_x$
& $-37 $ & $  18$
& $-36 $ & $  3$
\\
$\Delta F_x$
& $24 $ & $  5$
& $24 $ & $  1$
\\\hline\hline
\end{tabular}
\caption{Corrections to the chiral dispersive parameters in order to reproduce the KLOE data (in appropriate powers of $GeV$).}
\label{tab:cor}
\end{table}
Since we have no further information from KLOE concerning the efficiency
of every selected bin, the appointed errors in both tables are not very reliable.

As we have discussed above, it is not easy to interpret the result (moreover, if the fit relies just on the KLOE-like distribution we have made just from the 4+1 measured Dalitz parameters and do not have any notion of the systematic uncertainties here), but let us present some interesting observations. All the corrections of the parameters seems to reduce the original $O(p^6)$ contributions to these parameters coming from \cite{BG}, some of them even tend to change the sign of the total $O(p^6)$ contribution with respect to the original one. Taking into account relation (\ref{Eplus2F}) together with the condition (\ref{f rovna se g}) imply that in the case that all the difference between the physical amplitude and the original $O(p^6)$ chiral one is hidden in the wrong determination of the $C_i$s,
\begin{equation}
\Delta E_x+2\Delta F_x =0.
\end{equation}
From Table~\ref{tab:cor} we see this tendency.

Let us return to our main interest, the determination of the ration $R$. The integrated decay rate computed from this distribution is
\begin{equation}
 \Gamma(\eta\to\pi^+\pi^-\pi^0) = \sin^2 \epsilon \times 2.24(10) \text{ MeV}.
\end{equation}
Together with the experimental input for $\Gamma$ this leads to
\begin{equation}
 R  = 37.7(9).\quad
\text{[ChPT+disp.+KLOE]}
\end{equation}
The quoted error is based only on getting values for parameters in different regions as explained in the previous text.
It does not take into account systematic errors coming from the experimental data, which we do not know, and the theoretical error from the ChPT part. The later one can be estimated from the convergence of results coming from the chiral expansion.

The LO value coming from the current algebra is $R=19.1$ and the NLO result was $R=31.8$. Taking these values into account, our prediction for this ratio from the NNLO ChPT and data of KLOE is
\begin{equation}
 R  = 37.7\pm 2.8.\quad
\text{[ChPT+disp.+KLOE]}\label{Roo}
\end{equation}

\subsection{Overall fit: using the measured distribution}\label{sec:res}
So far we were using the dispersive representation in a very close connection with ChPT and up to the fact that we were trying to correct it by the values of the $C_i$s coming from KLOE, we just reproduced the amplitude
coming from this theoretical framework. But we can also change our strategy totally; we can use the experimentally measured distribution and employ ChPT just for the normalization. In this analysis we therefore assume that no matter what the proper description of this process leading to the correct physical amplitude would be, it would fulfill the general principles of quantum field theory together with the observed hierarchy of various contributions\footnote{This hierarchy is expressed in the construction of the representation in terms of the very basic chiral counting of the partial waves of the amplitude -- cf.~relations (2.2) and (2.3) in \cite{nas-proceeding2}.} which were used for the construction of our analytic dispersive representation. It then means that such correct physical amplitude can be to a good extent described by our parametrization and the later can be used for its analytic continuation to some region where ChPT gives a reliable result for the amplitude, and can be used there for the matching. Such procedure will substantially reduce the influence on $R$ of the error connected with the chiral expansion of the amplitude.

Let us start with the easier part --- fitting the KLOE-like distribution, which plays in this analysis a role of the experimental distribution. In contrast to our previous fits, where respecting the chiral orders of the dispersive parameters was natural (and important), in this case keeping the different chiral orders of the parameters makes no sense. The more natural approach is using our representation in the \emph{resummed} form --- the values of the parameters in the polynomial and in the unitarity part are the same.

The fit of this general representation to KLOE-like distribution was performed for the same data set as in the previous subsection (174 data points). Note that in the fit, the overall normalization is set so that the amplitude is equal to one at the center of the Dalitz plot. In order to simplify the comparison between this and the previous fits, we have decided to multiply all these data by the numerical factor that produces the same number for parameter $A_x$ as the corrected value obtained from the analysis in the previous subsection (average over all data sets), i.e.\ $A_x =
A_x^{(6)}+ \Delta A_x =  0.575$.

\begin{table}
\centering
\begin{tabular}{@{\hspace{1em}}c@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{2em}}r@{$\:\pm\:$}l@{\hspace{1em}}}\hline\hline
& \multicolumn{2}{c@{\hspace{2em}}}{cor.set} & \multicolumn{2}{c@{\hspace{1em}}}{fit to KLOE} \\\hline
$ A_x$
& $0.575 $ & $ 0.006$
& $0.575 $ & $ 0.001$
\\
$ B_x$
& $1.99 $ & $ 0.04$
& $2.15 $ & $ 0.02$
\\
$ C_x$
& $-6.8 $ & $ 0.3$
& $-5.8 $ & $ 0.2$
\\
$ D_x$
& $0.94 $ & $ 0.03$
& $0.87 $ & $ 0.08$
\\
$ E_x$
& $-31 $ & $ 3$
& $-19 $ & $ 9$
\\
$ F_x$
& $20 $ & $ 1$
& $21 $ & $  5$
\\\hline\hline
\end{tabular}
\caption{Dispersive parameters as free parameters fitted to KLOE (ultimate column; in
appropriate powers of $GeV$). For a comparison in the first column we have also summarized the corrected
parameters of the previous subsection.}
\label{tab:resummed}
\end{table}

The values of so normalized dispersive parameters coming from the overall fit of KLOE are presented in the ultimate column of Table~\ref{tab:resummed}. For the comparison there are displayed also the corrected values of these parameters from the previous subsection. Note that these two sets of values obtained from fit to KLOE correspond to two amplitudes with different unitarity parts (the first one contains the parameters of ChPT respecting their chiral orders, whereas the later contains in the unitarity part exactly those values of the parameters appearing in the polynomial part) and that in both cases the unitarity part forms an important part of the amplitude. The agreement between these two sets that is seen in this table is therefore quite interesting (the small disagreement is seen only for parameter $E_x$ which is, however, given with the biggest error in both approaches). Let us stress once more that the particular overall normalization was taken only to simplify this comparison and we still have to remember that values in Table~\ref{tab:resummed} are multiplied by, for the moment, unknown constant.

To set this overall normalization is, in our opinion, the main issue of the dispersive study for
$\eta\to3\pi$, so we try to be as cautious as possible. First of all, we will rely on a set of points rather than only on one point even though it would be sufficient for setting the normalization. Thus, we need to select the region of the points where we believe ChPT result. In order to achieve this task, we discuss the following articles which could be important for its selection:
\begin{itemize}
\item[\rim1)] Adler zero condition
\item[\rim2)] correspondence with order-by-order fit and convergence of chiral orders
\item[\rim3)] plateau argument
\end{itemize}
Let us explain them in detail.  The point \rim1) is connected with the $SU(2)\times SU(2)$ theorem
which was summarized in Sec.~\ref{sec:adler}. According to this
theorem the values of the charged amplitude at the points (\ref{adlerxi}), in
a small vicinity of $SU(2)\times SU(2)$ Adler zero, are $O(m_\pi^2)$ for $m_\pi\to 0$ and therefore protected
from being large. At $O(p^4)$ the zero of the real part of the
amplitude, namely, the point $s=u=1.4 m_{\pi}^2$ belonging to
the above set of protected points, possesses the following additional convenient
property. The slope of the amplitude (which is generally not
protected by that theorem) develops $O(p^4)$ corrections that are accidentally small. In the previous dispersive analysis \cite{Anisovich:1996tx}, this fact was the main motivation for matching of the dispersively constructed amplitudes exactly at this point. However, as we have learned from Sec.~\ref{sec:adler}, there is no guarantee that at the points near the Adler zero the chiral corrections to any calculated order are small.
Nevertheless, the points where the real part of the amplitude vanishes (which are often  called ``Adler zeros`` in this context too) can serve as good reference points, or benchmarks, of the individual chiral orders.
In such a way these points were used in the analysis of the NNLO ChPT calculations \cite{BG} with a result that the best convergence of their positions is observed on the line\footnote{Note that these points do
not belong to the set (\ref{adlerxi}), i.e.\ they are not close to the Adler zero in the strict sense of the $SU(2)\times SU(2)$ theorem.}
$t=u$. On the other hand, according to the same analysis, the point $s=u=1.4 m_{\pi}^2$ does not seem to be particularly stable with respect to the $O(p^6)$ corrections. The matching of the overall normalization at the
``Adler zeroes`` has however the advantage that at these points just the imaginary part of the amplitude is matched and therefore the uncertainty corresponding to the not well known $O(p^6)$ LECs is eliminated (or suppressed when we match in their vicinity).

The second article \rim2) can help us to reformulate the previous conclusion in different words.
The order-by-order fit should be by construction very similar to the chiral expansion. The
convergence of this expansion was crucial in \cite{Anisovich:1996tx} for setting
the matching point in $s=u$ around the $O(p^4)$ Adler zero ($s=1.4 m_\pi^2$). However, following the
detailed analysis of \cite{BG} we have to conclude that this choice becomes to be very
dangerous for matching at $O(p^6)$ (for $s=u$ the imaginary part at NNLO is even of the opposite sign when compared with NLO). Much better convergence when
coming from NLO to NNLO is seen for $t=u$, which we use in the following.

Finally \rim3) reflects the stability of the points within the given region or
cut. Studying some physical observable (for example the decay width) as a function
of the matching point one would expect a plateau behavior in the correct region.

Using the previous arguments we fit the normalization of subthreshold
parameters for $t=u$ cut (below the physical threshold) matching only the
imaginary part of ChPT amplitude where we interpolate between both variants,
the standard and the resummed one (cf.\ also Fig.~\ref{figtu} below).
Within the straightforward analysis one obtains
\begin{equation}
 \Gamma(\eta\to\pi^+\pi^-\pi^0) = \sin^2 \epsilon \times 2.25(40) \text{ MeV},
\end{equation}
which corresponds to the value
\begin{equation}\label{R disp}
 R  = 37.8 \pm 3.3.\quad
\text{[disp.+KLOE]}
\end{equation}

This number depends more strongly on the data than it was in the case of the order-by-order fit of the previous subsection. This is the reason why we have dropped ``ChPT'' in
its description even though one should remember that it enters the determination of this number through the normalization as explained above. The sources of errors are thus of the two types: the uncertainties connected with the experiment --- the uncertainties of the experimental data we fitted and the accuracy of their parametrization by our dispersive representation; and the uncertainties connected with the normalization procedure --- the error that is induced by the analytic continuation of the parametrization from the region where we have fitted the data to the region where we have matched with ChPT and finally, the error of the the determination of the values of the amplitudes in this region from ChPT, the error of the chiral expansion.

\begin{figure}
\begin{center}
\includegraphics{tuc.eps}
\end{center}
\caption{The real (blue) and the imaginary (red) part of the amplitude for $t=u$ line.
The dotted lines represent order-by-order fit (Set~3), the dashed ones stand for resummed fit constructed from it (with the values of all parameters equal to their $O(p^6)$ values from Set~3), and finally the solid lines reflect an overall fit corresponding to KLOE values, with normalization set to interpolate between the dotted and the dashed lines for the imaginary part. The vertical lines indicate the physical region.}
\label{figtu}
\end{figure}

From the way our parametrization was constructed and from the normalization procedure described above, it should be obvious that we have concentrated mainly on the reduction of the errors of the second type. Moreover, since we used only a distribution coming from 4+1 Dalitz parameters measured by KLOE without any information of the systematic uncertainties in different regions of the Dalitz plot, it should be obvious that these uncertainties prevail and the error quoted in (\ref{R disp}) corresponds solely to them. Provided we have more precise measurement of the distribution of the physical amplitude, we could quantify the uncertainty coming from fitting these data to our parametrization and a deeper analysis of the ones coming from the normalization procedure would be required.

To conclude the study on the charged decay $\eta\to \pi^+\pi^-\pi^0$ let us
summarize all the analyses performed in this and the previous
Sec.~\ref{sec:ObOC} and~\ref{sec:corOO} in one plot (Fig.~\ref{figtu})
focusing on the $t=u$ line.

\section{Neutral decay: \texorpdfstring{$\eta\to3\pi^0$}{eta -> 3 pi(0)}}\label{sec:neutral}
Before we present our results for the $\eta\to3\pi^0$ amplitude, let us shortly discuss the number of free parameters
of our dispersive parametrization for the $\eta\to3\pi$ decays that has to be determined from a fit.

The polynomial part of the charged amplitude was given in (\ref{polstu}). For the $\eta\to3\pi^0$ decay the situation is simpler since
{\allowdisplaybreaks[1]
\begin{equation}
  {\cal P}_0 (s,t,u) =  A_0 m_\eta^2 +
 C_0 \bigl[ (s-s^c)^2 + (t-s^c)^2 + (u-s^c)^2 \bigr]
 +   E_0
\bigl[ (s-s^c)^3 + (t-s^c)^3 + (u-s^c)^3 \bigr].
\end{equation}
}
These two decays are related by the 2-particle unitarity and so in the unitarity part of the $\eta\to\pi^+\pi^-\pi^0$ decay there appear 2 parameters from the neutral decay, whereas in the unitarity part of the $\eta\to3\pi^0$ decay there appear 4 parameters from the charged one. When one takes the full isospin breaking into account there is no further connection between these two decays and one needs to determine all these parameters appearing in the considered amplitude.
However, in the case we work in the leading order of the isospin breaking, relation (\ref{M0Mxiso}) bounds these two amplitudes together and all the parameters of the neutral decay can be expressed in terms of the charged parameters. The number of the parameters needed to be determined in the various studies in these two cases are given in Table~\ref{tab:pocty}.
\begin{table}
\centering
\begin{ruledtabular}
\begin{tabular}{lcc}
& $m_{\pi^\pm}=m_{\pi^0}$ & $m_{\pi^\pm}\neqm_{\pi^0}$
\\\hline
$\eta\to\pi^+\pi^-\pi^0$ & 6 & 8\\
$\eta\to3\pi^0$ & 5 & 7\\
$\eta\to3\pi$ & 6 & 9
\end{tabular}
\end{ruledtabular}
\caption{Numbers of free parameters of the dispersive parametrization that are needed to be determined in the various analysis. The first two lines describe the individual studies of the charged and the neutral $\eta\to3\pi$ decays, while the ultimate corresponds to the combined fit of both of them.}
\label{tab:pocty}
\end{table}

From the table and the present status of information we have from experiment on these amplitudes (as summarized in
Sec.~\ref{Dalitz plot parametrization}), it is obvious why we were focusing just on the charged decay modes of $\eta$.
In that case we have six unknown parameters in our dispersive formula (in the case $m_{\pi^\pm}=m_{\pi^0}$) which could be saturated by five known Dalitz parameters of this decay. On the other hand, although the neutral decay is theoretically much simpler (having less parameters and there is no P-wave contribution to the unitarity part), so far only one Dalitz parameter ($\alpha$) was measured for $\eta\to3\pi^0$. The procedure elaborated in the previous section will not be
thus very reliable in this case.

Working in the $m_{\pi^\pm}=m_{\pi^0}$ we can obtain the values of the neutral dispersive parameters from the values of the charged ones (using (\ref{M0Mxiso})), from which we can compute the neutral Dalitz parameters.

Doing so we get $\alpha = -0.042$ for the order-by-order correspondence (ChPT with an additional $O(p^6)$ polynomial) and $\alpha = -0.047$ for the overall fit (fit to KLOE with the normalization from ChPT). This two predictions can be put together to:
\begin{equation}
\alpha = -0.044 \pm 0.004.
\label{falpha}
\end{equation}
In these error bars the systematic uncertainties from KLOE are not included. This error also does not take into account the isospin corrections to the relation (\ref{M0Mxiso}).

The corresponding study of KLOE \cite{KLOE-etac} based on the charged decay mode leads to
the value $\alpha= -0.038(3)(12)$, which is also higher (in absolute value) than are the most precise values on this quantity (see Table~\ref{tabulka Dalitz eta0}, e.g. MAMI-C: $\alpha = -0.032(3)$). Note that the similar study done recently in~\cite{Schneider:2010hs} obtained the same prediction at the leading order in $\pi\pi$ rescattering taking using KLOE data: $\alpha^{(1)}_\text{NREFT} = -0.042(2)\tbinom{3}{5}$. Employing higher orders to this prediction leads to $\alpha_\text{NREFT}\approx -0.06$, the value even more off the current most precise experimental determination. However, these two values were really based on the KLOE data. The independent determination of \cite{Schneider:2010hs} for the parameter $\alpha$ is in much better agreement with the present most precise experimental value (see e.g.\ Table~\ref{tab:inputdalitz}).

Therefore, we have verified a deviation in the neutral Dalitz parameter $\alpha$ obtained from the measurement of the charged Dalitz parameters (using isospin relation (\ref{M0Mxiso})) and the one obtained from the direct
measurement of the neutral decay. This discrepancy can be attributed to the fact that we
do not have at disposal the exact KLOE distribution but only the simplified one. The second
reason is connected with the following disadvantage of the dispersive representation when connected to the isospin simplification $m_{\pi^\pm}=m_{\pi^0}$.
This limit forces us to use only one value for the two masses of pions. In the case of $\eta\to \pi^+\pi^-\pi^0$ we use an average value of the charged and the neutral one, whereas in the neutral case there is used naturally the neutral mass itself. In ChPT we are free to make this change without changing the other parameters of the amplitude (at least in principle LECs do not depend on the light quark masses).
This is, however, not true for the polynomial parameters used in the dispersive approach, which do depend on the masses.
Indeed, repeating the calculation of $\alpha$ done for the ``average'' pion mass we obtain the value of $\alpha$ that is reduced by approximately $10\%$ with respect to our result~(\ref{falpha}), i.e.\ it goes in a good direction to the independent measurements.

Similarly as in the previous section for the charged decay, we summarize our analyses for the neutral decay in one
plot (Fig.~\ref{figtuz}) again on the line $t=u$, now for absolute value squared of the amplitude.
\begin{figure}
\begin{center}
\includegraphics{tuz.eps}
\end{center}
\caption{The absolute value (squared) of the amplitude $\eta\to3\pi^0$ for $t=u$
line.
The notation is the same as in the previous graph: the dotted line represents order-by-order
fit (set 3), the dashed line denotes a resummed
fit (with $O(p^6)$ set~3 values), and the solid line stands for an overall fit
corresponding to KLOE value. With the vertical lines we have demarcated the physical region.}
\label{figtuz}
\end{figure}
Note the only slight change of the slope ($\alpha$) when passing from the order-by-order fit (or equivalently NNLO ChPT) to the resummed one (which is equivalent to a different resummation of the $O(p^8)$ chiral orders in ChPT). A bigger change in $\alpha$, connected even with a flip of the sign, is seen only after inclusion of the KLOE data.

Before closing our study of the neutral decay mode let us return back to the resummed correspondence introduced in Sec.~\ref{sec:res}. The major problem addressed there, the overall normalization, can be totally ignored for a
special quantity, the ratio of decay rates
\begin{equation}
r_\eta = \frac{\Gamma(\eta\to \pi^0\pi^0\pi^0)}{\Gamma(\eta\to \pi^+\pi^-\pi^0)},
\end{equation}
where such normalization simply drops out. Using the numbers from fit to KLOE in Table~\ref{tab:resummed} one obtains
\begin{equation}
r_\eta = 1.475\pm 0.015,
\end{equation}
the value that is in an excellent agreement with the most precise measurement of \cite{Achasov:2006dv}:
$r_\eta= 1.46(3)(9)$. We have verified that our result is stable against isospin corrections and is very close to the original ChPT value \cite{BG} $r_\eta^{NNLO} = 1.47$. Note that the prediction of NREFT \cite{Schneider:2010hs} is somewhat smaller $r_\eta = 1.40(1)(4)$.

\section{Conclusion}\label{sec:conclusions}

The main purpose of this paper was a demonstration of the possibilities of using the dispersive approach in the case of the $\eta \to 3 \pi$ decays. In the advent of the new precise measurements which could be sensitive to higher-order isospin-breaking effects it is a reliable theoretical tool, based only on general assumptions, as relativistic invariance,
unitarity, analyticity and crossing symmetry together with chiral counting, which ensures that our
amplitude is valid up to and including $O(p^6)$. It easily reproduces the ChPT amplitudes, providing them in a simple and compact analytic form. It is therefore useful for studies of the structure of the amplitudes and enables us to identify the basic independent combinations of the LECs appearing in them. In addition, it is more general and can be analytically continued outside the kinematic decay region. Further, it incorporates naturally the isospin corrections induced by the $m_{\pi^\pm}-m_{\pi^0}$ difference.

However, in the present situation when experiments are limited to 4+1 (one parameter is
zero or close to zero) Dalitz parameters for $\eta\to\pi^+\pi^-\pi^0$ and 1
parameter for $\eta\to3\pi^0$ decay it seems reasonable to use only leading
order of the isospin breaking. This leading order is hidden in the parameter $R$
(see (\ref{ARM})) and ${\cal M}(s,t,u)$ is thus effectively taken in the isospin
limit. Even such simplified analysis has an important theoretical outcome since by comparison of the experimentally measured decay rates of these processes with the predictions coming from theory, we can determine\footnote{Note that instead of $R$ we could use in our normalization of the amplitudes the parameter $Q$ as well. Such choice is favored in the analyses based on ChPT at NLO since at this order $Q$ depends only on QCD meson masses \cite{Gasser:1984gg} and is reasonably stable with respect to the Kaplan-Manohar transformation \cite{Kaplan:1986ru} of quark masses of ChPT. However, both these advantages of $Q$ are lost when one includes chiral two-loop effects because the relation between $Q$ and the meson masses gains noticeable $r=m_s/\hat{m}$ dependent chiral corrections at NNLO (cf.\ \cite{Amoros:2001cp}). When matching the amplitudes with the results of NNLO ChPT \cite{BG}, it is more natural to employ the normalization containing $R$ and assume that the Kaplan-Manohar ambiguity is fixed by the values of LECs used in that computation (for instance the value of $L_6^r=0$ stemming from large $N_c$ considerations).} the parameter $R$, thereby obtaining information on the individual masses of $m_u$ and $m_d$.

Such analysis was performed within NNLO ChPT, but the necessity of determination of the large amount of NNLO LECs of ChPT, together with the observed discrepancies between the values of the Dalitz parameters describing the energy dependence of the amplitude predicted by ChPT and those measured by experiments question also the accuracy of the result obtained for $R$.

On the analysis of the charged data from KLOE we have presented two methods how to make use both the information we have from the NNLO ChPT and the one from experiment together in order to determine the corrected value of this parameter. The first one is motivated by the possibility that all the discrepancy between the NNLO ChPT predictions and the KLOE measurement comes from the incorrect determination of the $O(p^6)$ LECs of ChPT (or more generally that all the discrepancy can be approximated well by a polynomial of the third order). Using it we have obtained the value $R = 37.7\pm 2.9$. The second method is based on the fit to the experimental data, in which the error of a slower convergence of the chiral amplitude is reduced by using ChPT just for the normalization of the parametrization in the region where the chiral expansion seems to converge well. Here the main source of uncertainties was due to the experimental error and we obtained $R=37.8\pm 3.3$. Even though both of these results depend on KLOE, the nature of the dominant error in each of the methods is different and so we can combine them to obtain our final prediction\footnote{Let us emphasize once more that this number stands and falls with the assumption that the genuine physical amplitude is described well by the KLOE-like distribution constructed from the values of 4+1 Dalitz plot parameters presented by KLOE \cite{KLOE-etac}. It is therefore desirable to remeasure the charged $\eta\to3\pi$ decay and repeat the performed analysis with the access to the real data.}
\begin{equation}
R=37.7\pm 2.2.
\end{equation}

This value agrees very well with the lattice average prediction~\cite{Colangelo:2010et} and is compatible with the Dashen's limit at NNLO order. (However, it would correspond to large Dashen violation if one worked at NLO as studied in \cite{Amoros:2001cp} --- cf.\ also Table 6 in \cite{BG}).
Using relation (\ref{RQ relace}) and the value of $r\equiv m_s/\hat m=27.4\pm0.4$ provided by lattice simulation~\cite{Colangelo:2010et}, we obtain
\begin{equation}
Q= 23.1\pm0.7.
\end{equation}
Note that \cite{Oller:2006xb} and several non-lattice determinations (for the recent overview see \cite{Leutwyler:2009jg,Zdrahal:2011ta}) point towards a smaller value for $r$, which will mean also a smaller value of $Q$, e.g.~the value of $r\sim24$ leads to $Q\sim21.5$. For the sake of completeness, let us also mention the value $Q=20.7\pm 1.2$ obtained from large deviation from Dashen's limit in \cite{Ananthanarayan:2004qk} (cf.\ also \cite{Kastner:2008ch}).

In the case we want to fully complete the task set in the Introduction and determine the current masses of $m_u$ and $m_d$ at some scale $\mu$, we need as inputs the values of $\hat{m}$ and $m_s$ at the same scale obtained from somewhere else. Introducing the most recent averaged values of PDG \cite{PDG}, in $\overline{\mathrm{MS}}$ scheme at the running scale $\mu=2\,\mathrm{GeV}$, $m_s=(100.6\pm2.1)\,\mathrm{MeV}$ and $\hat{m}=(3.8\pm0.1)\,\mathrm{MeV}$ into the definitions (\ref{hatm}) and (\ref{RQ}) as inputs, one obtains the current masses at this scale
\begin{equation}\begin{split}
m_u(\overline{\mathrm{MS}},\mu=2\,\mathrm{GeV})&=(2.52\pm0.13)\,\mathrm{MeV}, \\
m_d(\overline{\mathrm{MS}},\mu=2\,\mathrm{GeV})&=(5.08\pm0.13)\,\mathrm{MeV}.
\end{split}\end{equation}
Employing instead the averaged lattice results \cite{Colangelo:2010et} $m_s=(94\pm3)\,\mathrm{MeV}$ and $\hat{m}=(3.43\pm0.11)\,\mathrm{MeV}$ leads to
\begin{equation}\begin{split}
m_u(\overline{\mathrm{MS}},\mu&=2\,\mathrm{GeV})=(2.23\pm0.14)\,\mathrm{MeV}, \\
m_d(\overline{\mathrm{MS}},\mu&=2\,\mathrm{GeV})=(4.63\pm0.14)\,\mathrm{MeV}.
\end{split}\end{equation}
Our results are compatible both with the averaged values of these individual masses quoted in PDG \cite{PDG} and with the averaged lattice results of \cite{Colangelo:2010et}. Note that these two sets for each mass are correlated since the PDG values contains among others also the lattice results.

We have also made a quick look onto neutral decay mode $\eta\to3\pi^0$.
Assuming that the original KLOE measurement of charged Dalitz parameters
$a,b,d,f,g$ is correct we have verified some deviation in neutral Dalitz
parameter $\alpha$ obtained from our representation using them (and the isospin relation (\ref{M0Mxiso})) and from its direct measurements. In the previous section we have discussed the possible explanations of this discrepancy.
Note that the neutral decay $\eta\to3\pi^0$ hides more --- in the physical region there occurs an intrinsic cusp, however, its position is naturally very close to the edge of the phase-space. Its shape is
thus very suppressed and the techniques similar to the extraction of $K_L\to3\pi^0$
have to be employed. An unambiguous description of the amplitude including corrections
induced by $m_{\pi^\pm}-m_{\pi^0}$ difference can thus naturally help in this extraction.

At this point let us mention also the second neutral Dalitz parameter $\beta$, whose
measurement is planned also in the forthcoming experiments (either WASA or KLOE-II).
Its theoretical determination in ChPT is not influenced by $O(p^6)$ LECs and is
thus important for its consistency check. In Sec.~\ref{sec:Ci} we also list the other combinations of the Dalitz parameters that are suitable for such a check and their evaluation can shed light on the discrepancies in the determinations of the Dalitz parameters --- not only the discrepancy in $\alpha$ whose solution was proposed in \cite{Schneider:2010hs} but also the discrepancy in $b$ that even in \cite{Schneider:2010hs} remains.

The last physical quantity we have discussed was $r_\eta$, the ratio of the
neutral and charge decay width, with the result
\begin{equation}
r_\eta=1.475\pm 0.015.
\end{equation}
Our prediction, or more precisely the prediction of KLOE based on our dispersive parametrization, agrees not
only with the most recent experiment but also with the NNLO ChPT calculation.
Naively, seeing Figures \ref{figtu} and \ref{figtuz} and the change in the absolute value for both amplitudes, one would expect some impact, but apparently these two shifts canceled out in the ratio.

\acknowledgments{We would like to thank Hans Bijnens, Andrzej Kup\'s\'c and Bachir Moussallam for stimulating discussions and/or remarks on the text.
We are thankful to Hans Bijnens for providing us with \textsf{FORTRAN} programs which enable two-loop ChPT calculations.
Last but not least we have benefited from discussions with our late colleague Jan Stern.
This work is supported in part by the European Community-Research
Infrastructure Integrating Activity ``Study of Strongly Interacting
Matter'' (HadronPhysics2, Grant Agreement n.\ 227431) and the Center for
Particle Physics (project no.\ LC 527) of the Ministry of Education of the
Czech Republic.}

redpajama_set_name: RedPajamaArXiv