Extended chiral Khuri-Treiman formalism for $\eta\to 3\pi$ and the role of the $a_0(980)$, $f_0(980)$ resonances

Recent experiments on $\eta\to 3\pi$ decays have provided an extremely precise knowledge of the amplitudes across the Dalitz region which represent stringent constraints on theoretical descriptions. We reconsider an approach in which the low-energy chiral expansion is assumed to be optimally convergent in an unphysical region surrounding the Adler zero, and the amplitude in the physical region is uniquely deduced by an analyticity-based extrapolation using the Khuri-Treiman dispersive formalism. We present an extension of the usual formalism which implements the leading inelastic effects from the $K\bar{K}$ channel in the final-state $\pi\pi$ interaction as well as in the initial-state $\eta\pi$ interaction. The constructed amplitude has an enlarged region of validity and accounts in a realistic way for the influence of the two light scalar resonances $f_0(980)$ and $a_0(980)$ in the dispersive integrals. It is shown that the effect of these resonances in the low energy region of the $\eta \to 3\pi$ decay is not negligible, in particular for the $3\pi^0$ mode, and improves the description of the energy variation across the Dalitz plot. Some remarks are made on the scale dependence and the value of the double quark mass ratio $Q$.


Introduction
The physics of QCD in the soft regime is dominated by the phenomenon of spontaneous symmetry breaking because of the presence of three light quarks in the standard model. The low-energy dynamics can then be described accurately through an expansion built from a chiral effective theory (e.g. [1] for a recent review). This approach, which applies in both the Euclidean and Minkowski space-times is, to some extent, complementary to the purely numerical lattice simulation method. In the effective theory, however, part of the information on the non-perturbative QCD dynamics is contained as sets of values of the chiral coupling constants. These are not known, a priori, except for their order of magnitude [2] and must be determined as part of the probing of the effective theory.
In QCD, isospin breaking phenomena are driven by m d − m u , the mass difference of the two lightest quarks. For low-energy observables, an isospin breaking ratio which conveniently absorbs some of the next-to-leading (NLO) chiral coupling constants was introduced in ref. [3] In general, isospin breaking effects induced by electromagnetism are comparable in size to those proportional to m d − m u and their precise evaluation is made difficult by the poor knowledge of the associated chiral coupling constants [4]. In this respect, the η → 3π amplitude plays a special role because these electromagnetic contributions vanish in the SU(2) chiral limit [5] and are thus expected to be suppressed. This has been confirmed in the work of refs. [6,7] who evaluated the contributions proportional to e 2 m u , e 2 m d . A number of recent high-statistics experiments have studied the 3π 0 decay mode of the η [8][9][10][11][12] as well as the charged mode π + π − π 0 [13][14][15][16]. An extremely precise knowledge of the energy variation of the amplitudes squared across the Dalitz plot, which are traditionally represented by a set of polynomial parameters, has now become available. These accurate experimental results allow for stringent tests of the theoretical description of the amplitude which must obviously be passed before one attempts to determine Q. The Dalitz plot parameters derived directly from the NLO chiral amplitude, which was first computed in [17], are in clear disagreement with experiment. For instance, the prediction for the parameter α, involved in the neutral mode, has the wrong sign. The same problems, essentially, are found in the resummed expansion approach discussed recently in ref. [18]. The computation of the amplitude at the next-to-next-to leading (NNLO) chiral order was performed [19]. The comparison with experiment again fails if one assumes a simple naive model for the O(p 6 ) couplings C i (classified in ref. [20]) which are involved as six independent combinations. The η decay amplitude thus contains crucial information on the true QCD values of these couplings, which are essentially not known at present.
One obvious deficiency of the chiral expansion when calculating scattering or decay amplitudes in physical regions is the lack of exact unitarity (as emphasised e.g. in ref. [21]) which is restored gradually when going to higher orders. In the case of η → 3π, a large amount of work was devoted to the problem of estimating these unitarity, or final-state rescattering, higher order corrections [22][23][24][25][26][27][28].
In the present paper we reconsider, more specifically, the approach followed in refs. [22,24,25]. The main underlying assumption is that the chiral expansion of the η → π + π − π 0 decay amplitude should be optimally converging in an unphysical region of the Mandelstam plane in the neighbourhood of the Adler zero. The amplitude in the physical region is then deduced from a well-defined extrapolation procedure based on the analyticity properties of amplitudes in QCD, which utilise the set of dispersive equations derived initially by Khuri and Treiman [29] and perfected in the work of refs. [22,24,25]. These equations implement crossing-symmetry and unitarity in a more complete way than more naive loop-resummation approaches [26].
In previous work, ππ rescattering was assumed to be elastic. This is essentially exact in the physical η decay region. However, the dispersive formalism involves integrals over an energy range extending up to infinity. A property of the ππ scattering amplitude in the isoscalar S-wave is the sharp onset of KK inelasticity associated with the f 0 (980) scalar resonance [30,31]. Because of isospin violation, the ηπ → ππ amplitude actually exhibits a double resonance effect from both the f 0 (980) and the a 0 (980) scalars [32] near the KK threshold. Our aim is to propose a generalisation of the Khuri-Treiman formalism which takes into account KK inelasticity in the unitarity relations for both ππ scattering and ηπ scattering (which may be viewed as an initial-state interaction). Some approximations will be made, which simplify considerably the practical implementation, such that crossing-symmetry is maintained at the level of the η → 3π amplitudes but not in the amplitudes involving the KK channel. In this multichannel formalism, the double resonance effect is taken into account in the dispersive integrals and, furthermore, the construction of the η → 3π amplitude becomes valid in an extended energy region which includes not only the η decay region but also a portion of the ηπ → ππ scattering region. This will allow us to study how the energy dependence induced by the two 1 GeV scalar resonances propagate down to the low energy region and quantitatively affects the Dalitz plot parameters. The fact that these resonances could be influential at low energy was pointed out previously in ref. [33]. The plan of the paper is as follows. We first review in sec. 2 the derivation of the one-channel equations in which the amplitudes satisfy elastic ππ unitarity in the S and the P -waves. In sec. 3 we write the unitarity relations including the ηπ and the KK channels. A closed system of unitarity equations involves, besides η → 3π, isospin violating components of ηπ → KK, ππ → KK as well as KK → KK amplitudes. A multichannel set of Khuri-Treiman integral equations is defined such that the solution amplitudes satisfy these unitarity relations. We then discuss in sec. 4 the matching between the chiral expansion amplitudes and the dispersive ones. We adopt a simple approach which consists in imposing that the difference between the chiral NLO and dispersive amplitudes vanishes at order p 4 . This provides four equations which, in the single-channel case determine completely the dispersive amplitude provided one had introduced exactly four polynomial parameters in the Khuri-Treiman representation. This is generalised to the multichannel situation, in which one introduces 16 polynomial parameters. Finally, in sec. 5, the results on the Dalitz plot parameters are presented and some remarks are made on the determination of the quark mass ratio Q.

Khuri-Treiman equations in the elastic approximation
We remind below how the dispersion relations-based equations derived by Khuri and Treiman [29] for K → 3π decay can be generalised and applied to η → 3π, following [22,24,25].

Single-variable amplitudes for η → 3π
As initially demonstrated in the case of ππ → ππ (see [34]), amplitudes involving four pseudo-Goldstone bosons satisfy, in a certain kinematical range, an approximate representation in terms of functions of a single variable which have simple analyticity properties (see [35] for a review). In the case of η → 3π, and neglecting quadratic isospin breaking, there are three functions involved [24,25]: M 0 , M 1 , M 2 . We will follow the notation of ref. [25] and write the η → π + π − π 0 amplitude as with an overall factor 1 ǫ L which is proportional to the isospin breaking double quark mass ratio Q 2 given in eq. (1) The Mandelstam variables are defined, as usual, as and satisfy General analyticity properties imply that η → 3π decay and ηπ → ππ scattering are described by the same function in different regions of the Mandelstam plane. The corresponding physical regions are illustrated in fig. 1. The analogous representation for η → 3π 0 involves the two functions M 0 and M 2 only and reads The representations (2) and (6) are accurate in regions of the Mandelstam plane where the imaginary parts of the partial-wave amplitudes with angular momentum j ≥ 2 in the s, t or u channels are negligible (compared to those of the j = 0, 1 partial-waves). In the case of η → 3π or ηπ → ππ, this condition is satisfied in the range where s, t, u are sufficiently small compared with the masses squared of the tensor resonances, i.e. 1 It is convenient to formally factor out ǫ L but the amplitude is actually of the form: T η→3π = ǫ L A + ∆m 2 K B + e 2 C where e is the electric charge and ∆m 2 K is the physical K 0 − K + mass squared difference (see appendix A). |s|, |t|, |u| < ∼ 1 GeV 2 . This condition is also satisfied exactly by the amplitude obtained from the chiral expansion up to order p 6 [19]. This will prove very useful for writing matching conditions.
The functions M I (w) are analytic in w with a cut on the positive real axis: 4m 2 π ≤ w < ∞. Based on Regge theory, we expect that the functions M 0 (w) and M 2 (w) should not grow faster than w at infinity, while M 1 should be bounded by a constant. In the one-channel Khuri-Treiman framework, both M 0 and M 2 are usually assumed to behave linearly in w when w → ∞. With this asymptotic behaviour, there is a family of redefinitions of the functions M 1 and M 2 which leaves the physical amplitude A(s, t, u) unmodified [25,36], where z is the cosine of the scattering angle in the centre-of-mass frame of ηπ → ππ which is related to the Mandelstam variables by with From the representation of the isospin amplitudes, eq. (21) one easily derives the expression for the partial-waves. The result, for the j = 0, 1 partial-waves, can be written as, where the functionsM I (s) are given by linear combinations of angular integrals of the functions M IM with the notation [25] z n M I (s) = 1 2 Writing unitarity relations, one must first consider the unphysical situation where the η meson is stable, m η ≤ 3m π . The physical case is defined using analytic continuation in m η as in the classic derivation of generalised unitarity [37]. The contribution from the ππ states to the unitarity relation for the partial-wave M I j , reads with and f I j (s) is the ππ → ππ partial-wave amplitude, which is related to the scattering phase-shift by exp(2iδ I j (s)) = 1 + 2iσ π (s)f I j (s) .
The equality between the imaginary part and the discontinuity in eq. (28) holds when m η < 3m π . For the physical value of m η the right-hand side of eq. (28) continues to give the discontinuity across the unitarity cut, while the imaginary part must be deduced from the dispersive representation.

Khuri-Treiman equations in the elastic approximation
In the unphysical situation when m η < 3m π , the cuts of the functionsM I (w) are located in the region Re[w] < 4m 2 π such that eqs. (25) correspond to a splitting of the partial wave amplitudes into two functions which have a separated cut structure. When the η mass is increased to its physical value, the m 2 η + iǫ prescription must be used [37] and this insures that the complex cut of theM I functions, which approaches infinitesimally close to the unitarity cut in the region 4m 2 π ≤ s ≤ (m η − m π ) 2 , remains well separated from it (see fig. 4 in ref. [38]). Using the fact thatM I (s) has no discontinuity across the unitarity cut one can deduce from (28) that the discontinuities of the functions M I (s) along 4m 2 π ≤ s < ∞ are given by where j = 0 when I = 0, 2 and j = 1 when I = 1. In the sequel, we will drop the j subscript in the ππ phase-shift. Eqs. (31) imply that the functions M I can be written as Muskhelishvili-Omnès (MO) integral representations M 0 (s) = Ω 0 (s) α 0 + β 0 s + (γ 0 +Î 0 (s)) s 2 whereÎ , n a = δ 1a (33) and where the Omnès functions Ω I are given in terms of the ππ phase-shifts by the usual relation The phase-shifts, in the present context, are usually taken to obey the following asymptotic conditions δ 0 (∞) = δ 1 (∞) = π, δ 2 (∞) = 0 (35) which seem rather natural since these conditions are roughly satisfied at the KK threshold. In the elastic approximation framework, one can thus take the phases to be constant or quasi-constant in the inelastic region, above 1 GeV.
The polynomial part in the MO representation (32) was chosen to have four parameters. This allows one to define a unique dispersive amplitude by implementing four independent matching conditions with the chiral NLO amplitude. Taking into account the asymptotic conditions on the phase-shifts (35) one easily sees that the Khuri-Treiman equations (32) implement the following asymptotic behaviour for the functions M I The functions M I also satisfy the w = 0 conditions (10) and are thus uniquely defined.

Beyond the elastic ππ approximation
Elastic unitarity for ππ scattering is valid exactly below the four pions threshold and approximately up to the KK threshold. Close to 1 GeV, the ππ phase-shift increases very sharply under the influence of the f 0 (980) resonance which also couples strongly to KK. In order to properly account for the effect of this resonance it is thus necessary to go beyond the elastic unitarity approximation. We discuss in this section how to include both the ηπ and the KK channels into the unitarity relations and then generate a generalisation of the Khuri-Treiman equations. This will allow us to account for both the f 0 (980) and the a 0 (980) resonances in a realistic way.

ηπ contribution to unitarity
Including the ηπ channel in addition to ππ, the unitarity relation becomes where and f ηπ j (s) is the ηπ → ηπ partial-wave amplitude. In the energy region where ηπ scattering is elastic, which we will assume to extend up to the KK threshold, f ηπ j (s) is related to the scattering phase-shift by The j = 1 partial-wave f ηπ 1 corresponds to exotic quantum numbers j P C = 1 −+ and should thus remain rather small up to the 1 GeV region 2 . Therefore, ηπ rescattering is expected to affect mainly the two j = 0 amplitudes M 0 0 and M 2 0 . In the elastic regime, using eq. (37) one easily derives that the relation between the amplitudes on both sides of the unitarity cut reads, Comparing with the analogous relation in the elastic unitarity case (28) one deduces that including the effect of ηπ rescattering in the ηπ → ππ amplitude (which can be viewed as an initial-state interaction) amounts to simply perform the following replacements in the Omnès representations (32), In practice, ηπ rescattering is expected to become significant when the energy approaches the mass of the a 0 (980) resonance. It becomes necessary, then, to also take into account the KK channel.

KK contributions to unitarity
Let us now include the KK states into the partial-wave unitarity relations.
a) I = 1, j = 1: We are concerned mainly with j = 0 amplitudes but let us consider the j = 1 amplitude M 1 1 here also for completeness. The KK contribution reads The amplitude T ηπ + →K +K 0 1 is isospin violating 3 and, at linear order in isospin breaking, one can set σ K + K 0 (s) = σ K (s). We will denote the amplitudes appearing above as b) I = 2, j = 0: For the M 2 0 amplitude now, the KK contribution reads The amplitude T ηπ + →K + K 0 0 is isospin conserving in this case, since j is even while the amplitude T K +K 0 →π + π 0 0 is isospin violating. We will denote the two amplitudes in eq. (44) as 3 The amplitudes T ηπ + →K +K 0 j are isospin violating (conserving) for odd (even) values of j. This can be seen using G-parity: G|(KK) I j = (−1) I+j |(KK) I j . Since I = 1 for K +K 0 , G = +1 for odd values of j and −1 for even values while G = −1 for ηπ. c) I = 0, j = 0: Finally, for the M 0 0 amplitude, the KK contributions to the unitarity relations are Let us now separate the isospin conserving and the isospin violating contributions. For the kinematical factors, we introduce For the ηπ 0 → K + K − , K 0K 0 amplitudes, the isospin conserving part, g ηπ 0 , was introduced in eq. (45) and we call the isospin violating one N 0 0 . One has For the KK → (ππ) 0 amplitudes, the isospin conserving and the isospin violating amplitudes are denoted as Using this notation, the KK contributions in the unitarity relations for the partial-waves M I j can be summarised as follows These contributions involve new isospin-breaking KK → ππ and ηπ → KK amplitudes: In order to write a closed set of unitarity equations we must also consider KK → KK amplitudes. For these, the isospin conserving/violating components are denoted h I j , H 10 cons. viol.

Multichannel Khuri-Treiman equations 4.1 Closed system of unitarity equations
We can now write down a closed system of unitarity equations. It will be convenient to introduce a matrix notation for the isospin conserving amplitudes with I = 0 and I = 1 The I = 0 amplitude M 0 0 is now embedded into a system of four coupled unitarity equations where while the I = 2 amplitude M 2 0 is involved in a system of two unitarity equations, Finally, for the I = 1 amplitude M 1 1 , the coupled unitarity equations read In the following, however, we will disregard the inelasticity effects for M 1 1 and continue to use the elastic unitarity equation (28) in this case.

Coupled Muskhelishvili-Omnès representation
The next step is to write each one of the isospin violating partial-wave amplitudes as a sum of two functions, one having a right-hand cut only and one having a generalised left-hand cut. For physical mass values, the left and right-hand cuts of the various partial-waves appear to be overlapping, but it is possible to separate them unambiguously. In the case of the amplitudes involving the η meson, this is done by using the m 2 η + iǫ prescription. The amplitude KK → KK has a left-hand cut which extends on the real axis in the Using the m 2 K + iǫ prescription shifts this cut above the real axis.
For the I = 0 amplitudes one writes which generalises eq. (25) while for the I = 2 amplitudes one can write One can now employ the standard Omnès method in order to express the right-cut functions in terms of the left-cut ones. Introducing the matrix notation the discontinuity relation for the M 0 functions is deduced from the unitarity relation (53), where and ∆σ K is given in eq. (47). Eq. (60) generalises the one-channel discontinuity relation (31).
Let us now consider the following matrix where Ω I are the 2 × 2 Muskhelishvili-Omnès matrices corresponding to the T -matrices T I . Making use of the following discontinuity properties of the MO matrices one can express the discontinuity of the X matrix elements in terms of theM 0 functions where An alternative expression for ∆X a can be derived, using eqs. (63)) which shows that it represents the discontinuity of the following quantity across the right-hand cut. The quantity ∆X b is proportional to ∆σ K and it is given by We can then write a twice subtracted dispersive representation for X(s) and generate, via (62), a MO representation for the M 0 amplitudes where P 0 is a 2 × 2 matrix of polynomial functions, and the integral parts areÎ One remarks that in the term ∆X b the quark mass ratio ǫ L appears in the denominator and thus cancels with the overall factor in the complete amplitudes. This part is driven by the physical K + − K 0 mass difference via the ∆σ K function (see (47)). This mechanism was first studied in ref. [32] who predicted that a large isospin violation should take place at 1 GeV in the ηπ → ππ scattering amplitude, due to the contributions from both the a 0 (980) and the f 0 (980) resonances. The set of eqs. (69) account for the other sources of isospin violation as well.
We now consider the case of the I = 2 amplitudes and we also introduce a matrix notation for the column matrices From the unitarity relations including the ππ, πη and KK contributions we deduce the discontinuity of the M 2 functions As before, we introduce a matrix X 2 , multiplying M 2 by inverse MO functions, Its discontinuity relation is expressed in terms of theM 2 functions and the MO functions An alternative useful expression for ∆X 2 can be derived This leads to the following integral MO representation for the matrix M 2 , Eqs. (69) and (77) together with the uncoupled equation for M 1 (32) involve 16 polynomial parameters. The polynomial dependence was chosen such that the equations would reduce exactly to the set of elastic equations (32) if one switches off the coupling to the KK channel as well as ηπ rescattering 4 . One remarks that the asymptotic behaviour in the coupled-channel equations is modified as compared to the one-channel case: since the matrix elements of the MO matrices Ω I , decrease as 1/s when s goes to ∞ the entries of the M 0 matrix (thus the M 0 function) behave as constants. The asymptotic behaviour of M 2 , in contrast, remains the same as before.
In addition to the polynomial parameters, the right-hand side of eqs. (69), (77) involve a number of "hat functions". The functionsM I are determined in terms of the functions M I by the angular integrals (26), but one still needs to determine the other hat func-tionsĜ 10 ,N 0 ,Ĥ 10 ,Ĝ 12 which are related to the KK amplitudes. For this purpose, one would have to consider all the related crossed channel amplitudes an write similar sets of equations (which would introduce further one-variable functions). Here, since we are mainly interested in the ηπ → ππ amplitude we will content with an approximation for the amplitudes involving the KK channel, simply neglecting the integrals which involve the left-cut functions, i.e. we takê In support of this approximation, one observes that if one were to neglect the left-cut integrals in the η → 3π amplitude itself, one would still obtain a qualitatively reasonable description (e.g. [40]). With this approximation, eqs. (69) and (77)

Matching to the chiral amplitudes
We intend to fix the 16 polynomial parameters by matching to the chiral expansions of the amplitudes involved. For the η → 3π amplitude we will use the NLO expansion also including the part of the electromagnetic contributions of order e 2 (m d + m u ) from ref. [6]. This makes it possible to display explicitly the term induced by the K + − K 0 mass difference via unitarity which, in the dispersive representation, is contained in thê I b integrals (see eqs. (68), (71)). The explicit chiral expressions for the functions M I as used here are given in appendix A. For the isospin violating amplitudes involving the KK channel we will use the leading order chiral expansion. At this order, the partial-wave amplitudes have no left-hand cut, which is consistent with the approximation of dropping the left-cut functions in the integral equations. The relevant expressions (including the O(e 2 ) contributions) are given belowḠ and the coupling constant C can be related to the π + −π 0 mass difference (see appendix A). We implement matching conditions for the η → 3π amplitude following the simple method of ref. [38] which differs slightly from that of ref. [25]. LetM (s, t, u) be the amplitude computed in the chiral expansion at order p 4 . The polynomial parameters of the dispersive amplitude must be determined such that the dispersive and chiral amplitudes coincide for small values of the Mandelstam variables. At order p 4 one should thus have This condition is satisfied automatically for the discontinuities, which implies that one can neglect the discontinuity of the differences of the one-variable functions M I −M I and thus expand these differences as polynomials, (with n 0 = n 2 = 2, n 1 = 1). Inserting these expansions into the amplitude difference M(s, t, u) −M(s, t, u) and requiring that the O(p 4 ) polynomial part vanishes gives four equations. In the elastic Khuri-Treiman framework these determine the four parameters via the following four equations and one must keep in mind that the integralsÎ I (0) which appear on the right-hand sides depend linearly on the four parameters. In ref. [25] the first two of eqs. (84) were replaced by two equations related to the position s A of the Adler zero of the chiral amplitudē M (s, t, u) along the line u = s and the value of its derivative at s = s A . We will see below that eqs. (84) do actually implement these Adler zero conditions to a rather good approximation. Additionally, approximations were made in ref. [25] in the determination of β 1 and γ 0 (yielding e.g. γ 0 ≃ 0), the validity of which depends on the assumed behaviour of the phase δ 0 in the inelastic region. Here, the four equations will be solved without approximation. Doing so, one notes that the polynomial parameters get an imaginary part from the contributions of the integralsÎ I (0) which, however, is small (less that 10% of the real part).
In the coupled-channel case, the matching of the η → 3π amplitude again gives rise to four equations. In addition, we can match the values at w = 0 of each one of the four KK amplitudes with the chiral ones given in eqs. (80) as well as the values of the first and second derivatives (which are vanishing). Altogether, this provides 16 constraints which fix all the polynomial parameters appearing in the coupled-channel Khuri-Treiman equations. The details of these matching equations are given in appendix B.

Numerical method
The main difficulty which is involved in deriving accurate numerical solutions of the Khuri-Treiman equations is tied to the evaluation of the angular integrals z n M I needed to obtain the hat functionsM I and to the treatment of the singularities of these functions in the computation of theÎ a integrals which are finite. These technical aspects are discussed in detail in ref. [24]. By a change of variables one can rewrite the angular integrals as integrals over t with When s lies in the range (m η − m π ) 2 < s < (m η + m π ) 2 , the endpoints t ± (s) become complex. In fact, using the m 2 η + iǫ prescription one sees that t + and t − are placed on opposite sides of the unitarity cut when s gets larger that (m 2 η − m 2 π )/2. The integral from t − to t + must then be evaluated along a complex contour which circles around the unitarity cut of the functions M I . Rather than computing explicitly M I (w) for complex values of w, as in ref. [24], we follow here the approach of ref. [38] which consists in inserting the dispersive representations (11)  Convergence is found to be reasonably fast. Denoting by M . (88) Anticipating on the numerical results to be presented below, with n = 5, 6, 7 we find: ǫ (5) ≃ 4 · 10 −3 , ǫ (6) ≃ 2 · 10 −4 , ǫ (7) ≃ 4 · 10 −5 . The I = 2 amplitude is the one which has the slowest convergence rate.

Input I = 0, 1 T-matrices
Above the KK threshold, the S and T matrices are related by Two-channel unitarity implies that all the T -matrix elements are determined from three real inputs: a) the phase of S 11 , b) the modulus of T 12 and c) the phase of T 12 . For I = 0 scattering, experimental data exist for these quantities up to 2 GeV, approximately. We will use here a determination based on the experimental data of Hyams et al. [41] for the ππ → ππ phase-shifts and the data of refs. [42] and [43] (above 1.3 GeV) for the ππ → KK amplitude 5 . In the higher energy region, a smooth interpolation is performed which insure (in general) the existence of a corresponding unique MO matrix [45]. Below the KK threshold, a determination of the phase-shift based on the data of ref. [41] together with constraints from the ππ Roy equations (in the twice-subtracted version of ref. [46]) is performed. It is assumed that inelasticity can be neglected in this region, which implies that the phase of the ππ → KK amplitude coincides with the ππ phase-shift (Fermi-Watson theorem). This allows one to determine the modulus of this amplitude below the KK threshold [47] where it is also needed. Details on the parametrisation can be found in refs. [48,49]. Fig. 2 shows the experimental data used and the fitted curves.  In the case of the I = 1 T -matrix, experimental information on ηπ and KK scattering are indirect and far less detailed than for I = 0. We will rely here on the chiral Kmatrix model proposed in ref. [50] which provides a simple interpolation between certain known properties of the prominent I = 1 scalar resonances a 0 (980) and a 0 (1450) and the low energy properties constrained by chiral symmetry. The T -matrix reproduces, by construction, the ηπ → ηπ and ηπ → KK chiral amplitudes at NLO (and more approximately the amplitude KK → KK). It depends on the values of the O(p 4 ) chiral couplings: we use here (as in ref. [50]) a set of values for these taken from a p 4 fit of ref. [51]. We note that the value of L 3 in this set is L 3 = −3.82 · 10 −3 . We will use this value also in the computation of the η → 3π amplitude, for consistency.
A further constraint can be implemented by computing the ηπ and KK scalar formfactors from this T -matrix. Chiral symmetry relates the ηπ and the Kπ scalar radii, which leads to the prediction that r 2 ηπ S should be remarkably small, r 2 ηπ S ≃ 0.1 fm 2 . This small value can be reproduced provided the phase δ ηπ→KK raises sufficiently slowly above the KK threshold. The phenomenological parameters of the model are also constrained by the properties of the resonances. Fig. 3 shows a typical result for the ηπ phase shift and for the phase and modulus of the ηπ → KK amplitude 6 which we will use in the present work.
Finally, the phase shifts which we used for the I = 1 ππ P -wave and the I = 2 S-wave, for which inelasticity is ignored, are shown in fig. 4. In the energy region √ s ≤ 0.8 GeV, these phase shifts are given by the Roy equations solution parametrisation of ref. [46]. They are fitted to experimental data from [41] (P -wave) and [52,53] (I = 2) in the region √ s ≤ 1.5 GeV and interpolated to δ 1 1 (∞) = π, δ 2 0 (∞) = 0 in the higher energy region. c) In the lower energy region, s < ∼ 0.7 GeV 2 , on the contrary, the effect of the inelastic channels is to reduce the size of the amplitude. One also observes that in this region the influence of the inelastic channels becomes small. d) In the sub-threshold region, finally, the coupled-channel and single-channel amplitudes are essentially indistinguishable. This is expected since the amplitudes are constrained to satisfy the same chiral matching equations. 6 The phase shown is that of T 0 (ηπ + → K +K 0 ) = −g ηπ 0 (according to the isospin convention of eq. (13)). It satisfies δ 0 ηπ→KK (m a0(1450) ) = 100 • which corresponds to r 2 ηπ S = 0.12 fm 2 , within 20% of the chiral O(p 4 ) result. It is not difficult to identify the main mechanism which generates the behaviour described above. For this purpose, let us consider the i = 1, j = 1 component of the matrix M 0 and let us absorb the the effect of the integralsÎ a ,Î a at s = 0 into the polynomial matrix, definingP

Illustration of the role of the inelastic channels
The real parts of the three components ofP 0 which are related to the KK channel have the following expressions  These polynomial coefficients are controlled, we recall, by the matching conditions to the LO chiral ηπ → KK, ππ → KK, and KK → KK isospin violating amplitudes (see appendix B). The component Re [P 0 ] 21 is negative in the region s > 4m 2 π and dominates the others in the range s > ∼ 0.2 GeV 2 . In fact, the corresponding contribution to M 0 dominates in the whole region 4m 2 π < s < 0.95 GeV 2 . In this region, the contribution from the inelastic channels can thus be written approximately as The components of the Omnès matrices which appear in eq. (93) are plotted in fig. 6. The salient feature here is that the real part of the I = 0 component (Ω 0 ) 12 is positive at low energy and changes sign 7 at s ≃ 0.73 GeV 2 . This is the main reason why the inelastic channels decrease the η → 3π amplitude below 0.7 GeV 2 and increase it above. This behaviour is enhanced by the I = 1 component (Ω 1 ) 11 which is larger than 1 (reflecting the attractive nature of the ηπ → ηπ interaction below 1 GeV). The behaviour of the amplitude along the t = s (or u = s) line, which displays the Adler zero, is shown in fig. 7. In the sub-threshold region, the chiral, the single-channel and the coupled-channel amplitudes are seen to be very close. For the position of the Adler zero, we find Finally, the results for the η → 3π 0 amplitude are shown in fig. 8. The influence of the inelastic channels are seen to be quite substantial in this case in the whole low-energy region. One also sees that there is no region in which there is close agreement between the dispersive and the chiral O(p 4 ) amplitudes. This, of course, is because of the occurrence of the combination M 0 (s) + M 0 (t) + M 0 (u) and the fact that at least one of the variables 7 The presence of a zero below the KK threshold follows from Watson's theorem which is obeyed by the component (Ω 0 ) 12 and leads to the relation Re (Ω 0 ) 12 /Im (Ω 0 ) 12 = cot δ 0 0 . This implies that Re (Ω 0 ) 12 vanishes when the phase shift δ 0 0 goes through π/2. s, t, u must lie above the ππ threshold, thus generating significant S-wave rescattering chiral corrections.

Dalitz plot parameters
One traditionally describes the Dalitz plot in terms of two dimensionless variables X, Y such that |X|, |Y | ≤ 1 and the centre of the Dalitz plot corresponds to X = Y = 0. In the case of the charged decay amplitude η → π + π − π 0 , the variables X, Y are related to the Mandelstam ones by with Q c = m η − 2m π + − m π 0 . Assuming charge conjugation invariance, the amplitude must be invariant under the transformation X → −X and a polynomial parametrisation of the amplitude squared can be written as In the case of the neutral decay amplitude η → 3π 0 , Q c must be replaced by Q n = m η − 3m π 0 in the definition of X and Y . Charge conjugation and Bose symmetry imply that the amplitude must be invariant under the two transformations with z = X + iY . The amplitude squared can thus be represented as A direct comparison of the dispersive amplitudes squared with the experimental data from KLOE [16] is shown in fig. 9 and our numerical results for the Dalitz plot parameters are collected in table 2. The numbers quoted in the table are obtained in a way which parallels the experimental determination: a discrete binning of the Dalitz plot is performed (with a few hundred bins) and a global least squares fit of the theoretical (chiral or dispersive) amplitudes squared is performed using the representations (96), (98). The table also shows the two most recent experimental determinations [15,16].  Table 2: Results for the Dalitz plot parameters of the charged and neutral η → 3π decays based on the NLO chiral amplitude (in the form given in appendix A) and its dispersive extrapolations based on single-channel and coupled-channel Khuri-Treiman equations and a matching procedure described in the text. The last two columns show experimental results from refs. [15,16] and [54].
It is clear, at first, that the amplitudes obtained from solving the Khuri-Treiman equations and constrained to match the chiral NLO amplitudes are in much better agreement with the experimental results in the physical decay region than the NLO amplitude itself. In particular, the parameter b which was too large by a factor of three is reduced by a factor of two and the parameter α which was positive becomes negative. This is in close agreement with the results obtained a long time ago by Kambor et al. [24]. Our main new result is that taking into account the KK inelastic channels and the effect of ηπ rescattering has a non negligible influence on the Dalitz parameters and tends to further improve the agreement with experiment. The influence of these inelastic channels for the parameters d and g is small (less than 5%) but quite significant for the parameter b which is reduced by 17% and now lies within 15% of the experimental value. This reflects the reduction of the amplitude caused by the inelastic effects at low energy discussed in sec. 5.3. Similarly, the parameter α is modified by approximately 20% by the KK inelastic channels and becomes rather close to the experimental result. The parameter g is the only one which shows a mismatch, by a factor of two, with the value measured by KLOE.
Finally, let us consider the sensitivity of the Dalitz plot parameters to the strength of the ηπ interaction, which is not precisely known at present. This is illustrated in table 3. The table also shows that varying the O(p 4 ) coupling L 3 has a significant influence, in particular on the parameter d.

5.5
The ratio Γ(η → 3π 0 )/Γ(η → π + π − π 0 ) Let us quote here the results for the ratio of the 3π 0 and the π + π − π 0 decay rates We find that the influence of the inelastic channels on this ratio is very small, As compared with the chiral O(p 4 ) result this ratio is thus predicted to increase under the effect of the final-state interactions, by only 2%. The experimental status of this quantity is not completely clear at present, as the PDG [54] quotes two different numbers Besides, the CLEO collaboration [56] has performed an experiment dedicated to the determinations of the η meson decay branching fractions and they quote as the most precise determination of the 3π 0 /π + π − π 0 ratio.

The quark mass ratio Q from the chirally matched dispersive amplitude
It must first be reminded that, once the electromagnetic interaction is taken into account, the quark mass ratio Q is no longer invariant under the QCD renormalisation group since the quark mass variation with the scale depends on its electric charge (with e u = 2/3, e d = e s = −1/3). This implies the following scale variation for the factor It can then easily be verified, using the equations of appendix A which include the e 2 m 2 π contributions [6], that the scale invariance of the complete NLO chiral amplitude is restored thanks to the combination of two of the electromagnetic coupling constants [4], K r 9 + K r 10 . Indeed, as shown in refs. [57,58] this combination depends not only on the chiral scale µ but also on the QCD scale µ 0 and satisfies In practice, this means that in order to determine Q(µ 0 ) from the η → 3π amplitude we must specify the values of the electromagnetic chiral couplings K i at the corresponding scale. We choose here µ 0 = 0.77 GeV and estimate the values of the couplings K r i (µ, µ 0 ) from a resonance saturation model [59]. Such estimates are qualitative at best but it can be shown, based on general order of magnitude arguments, that the uncertainty induced on the amplitude should not exceed a few percent [6,7]. Having verified that the dispersive amplitude is in qualitative agreement with experiment concerning the Dalitz plot parameters, we can make an estimate of the value of Q. In the present approach, the amplitude in the physical decay region is derived as a Khuri-Treiman solution uniquely defined from the chiral NLO amplitude by the set of four matching equations, and thus has a definite dependence on Q. We can then estimate the value of Q by the requirement that the dispersive amplitude reproduces the experimental values of the η → 3π decay widths taken from the PDG [54] (constrained fit). Doing this, we find η → π + π − π 0 : Q = 21.8 ± 0.2 (single-channel), Q = 21.6 ± 0.2 (coupled-channel) η → 3π 0 : Q = 21.9 ± 0.2 (single-channel), Q = 21.7 ± 0.2 (coupled-channel) (108) where the quoted errors only reflect the experimental ones on the widths (we comment below on the theoretical error). This shows that the effect of the inelastic channels on the determination of Q is rather small, of the order of 1%, and tends to decrease its value.
The central value of Q is somewhat smaller than the results which are obtained from lattice QCD+QED simulations of hadron masses: Q = 22.9 ± 0.4 (ref. [60]), Q = 23.4(0.4)(0.3)(0.4) (ref. [61]). The error on Q associated with the phase-shifts below 1 GeV is rather small, of the order of 1%. The error associated with the NLO amplitude, essentially related to L 3 , is of the same order. The largest error arises from chiral NNLO contributions to the amplitude which will modify the determination of the polynomial parameters via the matching conditions. Assuming that they induce a 10% relative error in the determination of each one of the four polynomial parameters α 0 , β 0 , γ 0 , β 1 and assuming the errors to be independent, gives the following theory error on Q Q being mostly sensitive to the variation of the first two parameters α 0 , β 0 . We have also varied the remaining 14 polynomial parameters, assuming a 20% relative error, and found that they have a much smaller influence. Within the error (109), the determination based on η → 3π decay is compatible with the lattice QCD results, which confirms that the size of the NNLO corrections to the η decay amplitude in the sub-threshold region should not exceed 10%.

Further experimental constraints on Q
Our estimate for the error on Q (eq. (109)) was based on a general order of magnitude assumption on the size of the NNLO corrections to the four leading polynomial parameters. More precise information on the size of these corrections can be derived by making use of the precise experimental results on the energy dependence across the Dalitz region, imposing that the dispersive amplitude reproduces these via a least-squares fit. Not all the four leading polynomial parameters can get independently constrained in this way since a ratio of amplitudes is involved. We will make the simple choice to fix one of them, α 0 , from its chiral value and to perform a variation of the three others β 0 , γ 0 , β 1 . We will use the latest KLOE experimental data [16] which consist of a set of amplitudes squared, |T exp (X i , Y i )| 2 , measured over 371 energy bins in the Dalitz region and satisfying the normalisation condition We introduce corresponding theory amplitudes T th (X, Y ) = M c (X, Y )/M c (0, 0.05) and define the χ 2 as allowing for a floating of the normalisation within the experimental error via a parameter λ.
At first, setting λ = 1 and computing the χ 2 with our chirally matched central amplitude with L 3 = −3.82 · 10 −3 we obtain χ 2 = 3079. If, instead, we use the value recently derived in ref. [55] from K l4 decays: L 3 = −2.65 · 10 −3 , one obtains a significantly reduced result: χ 2 = 714. It thus seems reasonable to use this value as a starting point, which fixes the central value of α 0 . We have then searched for a minimum of the χ 2 by varying the real parts of the three polynomial parameters β 0 , γ 0 , β 1 (keeping their imaginary parts fixed) and the normalisation λ. The other polynomial parameters, in particular α 0 , are kept fixed to their chirally-matched value. In this way, we obtain at the minimum: χ 2 min = 387, which corresponds to a value per degree of freedom χ 2 min /dof = 1.055. The fitted values of the polynomial parameters differ from the chirally-matched values by less than 5% (the numerical values are given in appendix D). This confirms that this difference can be consistently attributed to chiral NNLO effects.
The theory error in this approach is dominated by the variation (by 10%) of the single polynomial parameter α 0 (we have also added an error associated with the input T -matrices and an error from varying the imaginary parts of the polynomial parameters). The error induced by the variation of the parameters β 0 , γ 0 , β 1 is now computed using the covariance matrix as evaluated by the MINUIT fitting program [62] (this matrix is given in appendix D). This error is added in quadrature with the experimental error on the π + π − π 0 and 3π 0 decay rates. Finally, the result for the quark mass ratio as determined from this fitted amplitude can be written as We find it useful to quote the theoretical and experimental errors separately since the former is not necessarily gaussian. Previous determinations of Q which combine chiral constraints and fits to Dalitz plot data have been performed in refs. [28,[63][64][65][66][67]. In the most sophisticated of these approaches [67], five polynomial parameters are included in the fit and the effect of the π + − π 0 mass difference in the amplitudes is accounted for. Finally, we quote the values of the 3π 0 Dalitz plot parameters which can be predicted from our fitted amplitude

Conclusions
We have proposed an extension of the Khuri-Treiman formalism for the η → 3π amplitude which includes the ηπ and the KK channels in the unitarity equations in addition to the elastic ππ channel. Modulo some approximations (in particular we do not attempt to impose unitarity in the crossed channels involving kaons like πK → πK or ηK → πK) the equations for the one-variable functions M 0 and M 2 are shown to be simply replaced by 2 × 2 matrix equations. These are given in eqs. (69) and (77) which involve both the I = 0 and the I = 1 Omnès 2 × 2 matrices. Eq. (69) exhibits explicitly the contribution induced by the physical K 0 − K + mass difference via unitarity in integral form. The amplitudes derived from this extended framework should be valid in an energy range which covers the physical decay region and also the physical region of the scattering ηπ → ππ below 1 GeV. Given a fixed number of polynomial parameters, an improved precision at low energy should result from the fact that the effects of the two prominent light scalar resonances a 0 (980) and f 0 (980) are taken into account in the dispersive integrals.
Using four polynomial parameters in the η → 3π amplitude we have reconsidered the idea of performing a prediction of the amplitude in the physical region as an extrapolation of the O(p 4 ) chiral amplitude, uniquely defined by fixing all the polynomial parameters by matching conditions. These are imposed in the form of a set of equations which insure that the differences between the dispersive and the O(p 4 ) chiral η → 3π amplitude are of order p 6 or higher. One verifies then, that the chiral and the dispersive η → π + π − π 0 amplitudes are very close in the neighbourhood of the Adler zero. These conditions also insure, for the charged decay amplitude, that the single and multi-channel dispersive amplitudes are quasi-identical in the whole region 0 ≤ s ≤ 4m 2 π , |t − u| ≤ (m η + m π ) 2 . In contrast, for η → 3π 0 , one finds that the unitarity induced chiral corrections are significant even in the sub-threshold region.
We have considered the Dalitz plot parameters and we found that the induced influence of the a 0 (980), f 0 (980) resonances is not negligible, in particular for the neutral mode. The modifications of the parameters, in the coupled-channel framework, go in the sense of improving the agreement with experiment, in particular for the parameters a, b, f of the charged mode. The parameter α, for the neutral mode is modified by 20% by the effects of the resonances and lies rather close to the experimental value. The remaining differences between the experimental and the dispersive-theoretical amplitude suggest that NNLO contributions are needed in the matching conditions, at the 5-10% level, which seems quite plausible. Some of these NNLO effects could be accounted for in a more general framework which would implement both unitarity and crossing symmetry completely for all the channels involved. This is left for future work.
The η → 3π amplitudes constructed in the present approach inherit a well defined dependence on the quark mass ratio Q from that of the chiral NLO amplitude. We can then determine Q such as to reproduce the integrated decay widths. The central value that one obtains is somewhat low compared to the recent determinations from lattice QCD simulations but it is compatible within the uncertainty induced by the NNLO effects in the matching. Some knowledge of these NNLO effects seems necessary in order to improve the precision of the determination of Q by this method. two physical quantities 8 and in terms of the quark mass ratio (m s +m)/m (see ref. [3]). We also include the electromagnetic contributions of order e 2 p 2 evaluated in ref. [6] which allows one to identify the piece induced by the physical K + − K 0 mass difference via unitarity. Further electromagnetic corrections which have been computed in ref. [7] are not included here. The expressions forM I given below also implement the w = 0 conditions (10), which simplifies somewhat the writing of the matching relations. Using the following notation the functionM 0 reads: The contributions ∆M a 0 , ∆M b 0 are induced by the electromagnetic interaction and proportional to e 2 , The last contribution, ∆M c 0 , is induced by the physical The parameters C and K r i in the above expressions are the chiral coupling constants which appear at order e 2 and e 2 p 2 respectively [4].

B Matching equations
We reproduce below the set of matching relations from which we determine the set of 16 polynomial parameters of the Khuri-Treiman coupled-channel equations (i.e. eqs. (69), (77) and the second one of eqs. (32)). In order to simplify the relations it is assumed here that the chiral expressions for the η → 3π functionsM I satisfy, as in appendix A, the w = 0 relationsM 1 (0) =M 2 (0) =M ′ 2 (0) = 0. Derivatives at w = 0 are denoted either by dots or by primes and matrix elements of the I = 0, 1 MO matrices are denoted here by Ω

C Angular integrals and hat functions
Using the dispersive representations (11) of the functions M I , one can express the angular integrals in the following form which displays explicitly their singularity when s → (m η − m π ) 2 . For I = 0, 2 one has M I (s) =α I + 1 2β where κ(s) is given in eq. (23). For I = 1 one has The functions R n I (s) which control the singularities arise from the part of the t integration contour which encircle the unitarity cut, they are given by (where t ± (s) is given in eq. (87)) in the s range and R n I (s) = 0 otherwise. In particular, no divergence occurs when s → 4m 2 π or s → (m η + m π ) 2 .