Three-body final state interaction in $\eta \to 3 \pi$ updated

In view of the recent high-statistic KLOE data for the $\eta \to \pi^+ \pi^- \pi^0$ decay, a new determination of the quark mass double ratio has been done. Our approach relies on a dispersive model that takes into account rescattering effects between three pions via subenergy unitarity. The latter is essential to reproduce the Dalitz plot distribution. A simultaneous description of the KLOE and WASA-at-COSY data is achieved in terms of just two real parameters. From a global fit, we determine $Q=21.6 \pm 1.1$. The predicted slope parameter for the neutral channel $\alpha=-0.025\pm 0.004$ is in reasonable agreement with the PDG average value.


I. INTRODUCTION
Three meson systems play an important role in studies of hadron reaction dynamics and hadron spectroscopy. For example, in three-particle decays of heavy quarkonia several candidates for non-quark model resonances have recently been observed [1][2][3]. Three-body decays of B and D mesons are a promising laboratory for studies of CP-violation [4,5]. In the light meson sector the limited phase space makes three-particle decays an ideal testing ground of effective theories of strong interactions. Detailed amplitude analysis of three meson production becomes even more important in light of the current and forthcoming high precision data from various hadron facilities [6][7][8][9].
The isospin breaking η → 3π decay, which we consider here, is of great importance as it allows to measure the light quark mass difference. The electromagnetic effects are known to be small [10][11][12], and the decay is driven by strong interactions through the ∆I = 1 isospin breaking transition that appears directly in the QCD Lagrangian. The decay amplitude is proportional to the light quark mass difference, (m u − m d ), and it is conventionally expressed in terms of the parameter Q 2 defined by with m s being the strange quark mass. Note, that this double ratio (1) is protected to strong high order corrections. Given the small breakup momenta the distribution * Electronic address: danilkin@uni-mainz.de of pions in the Dalitz plot of the η → π + π − π 0 decay can be analyzed in terms of a small number of parameters that determine deviations from a uniform distribution. These parameters are referred to as Dalitz plot parameters. Early analyses [13][14][15] could determine only a few, leading Dalitz plot paramters. A few more parameters were determined using the 2008 KLOE measurement [16]. These analyses were further improved thanks to the highquality WASA-at-COSY [17] and new KLOE [18] data. The statistics of the most resent measurements is high enough to allow for binned, data-driven analysis. On the theoretical side there has been significant progress in chiral, effective field theory analysis of η decays. Chiral perturbation theory (χPT) seems to converge poorly, yielding Γ η→π + π − π 0 = 66, 167 ± 50, ∼ 300 eV at leading (LO), next-to-leading (NLO) and next-tonext-to-leading (NNLO) order, respectively [19][20][21][22]. This indicates importance of pion-pion interactions and various approaches have been used to implement these effects to all orders [23][24][25][26][27].
In our recent study [28] we implemented the S-matrix constraints of unitarity, analyticity and crossing symmetry via a set of dispersion relations [29][30][31]. Connection with QCD is achieved by matching the dispersive amplitudes with χPT at the point where the latter converges best i.e. below the threshold. In [28] we used WASAat-COSY data [17] to determine the free parameters, i.e. subtraction constants of the dispersive integrals. As a result, we achieved a simultaneous description of the Dalitz plot distributions of the charged and neutral η decay modes. In order to extract the parameter Q we matched the dispersive amplitude with the next-to-leading order (NLO) χPT result near the Adler zero and we obtained Q = 21.4±0.4 [28]. The purpose of this letter is to revisit the result of [28] in a view of the new high statistic data from the KLOE experiment [18].

II. THE METHOD
In this section, we briefly review the η decay amplitudes that were developed in [28].
For η → 3π, the transition amplitude A(s, t, u) is a function of three Mandelstam variables s = (p π + + p π − ) 2 , t = (p π − + p π 0 ) 2 , and u = (p π + + p π 0 ) 2 which are related by s + t + u = m 2 η + 3 m 2 π . Except for the phasespace boundary, we work in the isospin limit and take m π = (2 m π + + m π 0 )/3. At low energies, one can perform a partial wave (p.w.) decomposition while crossing symmetry implies unitarity cuts in all three Mandelstam variables. Therefore we symmetrize the p.w. expansion in all three channels [29,[32][33][34][35][36][37], which for the charged decay, η → π + π − π 0 , implies the following representation, The amplitudes a IL have only the right-hand, unitary cuts. In Eq. (2) z i ≡ cos θ i and θ s,t,u are the centerof-mass scattering angles in the s, t and u-channels, respectively. The subscript (I, L) labels isospin and orbital angular momentum, with I + L = even due to Bose symmetry. The latter implies that for L max = 1 there are three unknown isospin amplitudes with (I, L) = (0, 0), (2, 0), (1,1). An amplitude of the neutral decay, η → 3π 0 can be easily reconstructed using ∆I = 1 relation, We emphasize, that the decomposition in Eq.(2) has the same analytical properties as the amplitude in NNLO chiral expansion [38,39]. However, in contrast to χPT, we can impose unitarity to all orders on the a IL amplitudes. The discontinuity along the right-hand cut can be expressed through the elastic ππ partial wave amplitudes f IL , which leads to Note, that the amplitudes a IL (s) and f IL (s) for L > 0 are subject to kinematical constraints [40] which have to be removed before application of dispersion relations. This is done by introducing the reduced amplitudes a IL (s) = a IL (s)/Z L (s) where the factor Z L (s) is proportional to the product of the c.m. momenta of ππ and πη [28]. The normalization is fixed by the phase space factor ρ(s) = 1 − 4 m 2 π /s, with Im(1/f IL (s)) = −ρ(s). The explicit form of the crossing matrices C II st,su can be found in [28]. The contribution from the first term on the righthand side of Eq. (4) reproduces the direct s-channel unitarity, while the second term contains the left-hand cuts from the t and u-channels. While calculating the latter, special care has to be taken for 4 m 2 π ≤ s < (m η + m π ) 2 , i.e. one has to deform the contour to avoid the cut along the real axis [32,41,42]. The kinematical singularity free amplitudesã IL (s) satisfy the Cauchy representation, which up to subtraction constants yields, The combination of (5) and (4) sets the so-called Khuri-Treiman (KT) framework that can be solved using several techniques [43]. The most popular method is to write a set of dispersion relation for the ratioã IL (s)/Ω IL (s) with Ω IL (s) being the Omnès function [23,27]. In this case it is necessary to make further assumptions about the unknown high-energy region which is typically done by introducing subtractions. This procedure is relatively easy for ω → 3π decay which depends dominantly on the pion-pion P-wave scattering input [44,45]. However, this is more challenging for η → 3π decay, where the dominant contribution comes from the S-wave and the form of the ππ isoscalar Omnès function is very sensitive to the asymptotic behavior of the phase shift δ IL=00 (s → ∞).
In [46] (Figs. 4 and 7) different scenarios for the the ππ phase shift inputs were investigated and the corresponding Omnès functions were produced. In order to minimize these differences, the subtraction polynomial of the sufficient order is required in the dispersion representation. A complementary approach is the Pasquier inversion [36,47] that we applied to analyze WASA-at-COSY data in [28] and use here as well. This method uses contour deformation to exchange the order of double integral appearing on the right hand side of Eq. (5). As it was shown in [43], once the two-body amplitudes f IL (s) are known, different methods for solving the disersive integral provide the same result. However, when the Pasquier inversion is applied, the input of f IL (s) is required in a different energy region.
To proceed we writeã IL (s) in the form where the function F IL (s) is introduced to remove Adler zeros specific to the elastic amplitude f IL (s) and introduce zeros in the decay amplitude as required by chiral symmetry. From Eq. (4) and Eq. (6) one can derive the discontinuity of g IL and write the dispersive representation for g IL (s). As a result we obtain a double integral equations for g IL (s), which can be reduced to a single integral equation using the Pasquier inversion, The explicit form of the kernel functions, K IL,I L (s, t) can be found in [28]. Currently, they are only calculated in the region s ∈ (0, (M − m π ) 2 ), which includes physical region and therefore cover the whole Dalitz plot region. In order to compute the amplitudes beyond that region, a proper analytical continuation is required within Pasquier inversion technique. Important steps in that direction were already elaborated in [48]. We will come this issue later in the Q-value determination.
The first term and the part of the second term on the right-hand side have the left-hand cut and can be expanded in the Taylor series in the physical region. Retaining only a single term in the expansion we arrive at the following relation [43,49] which is solved by discretizing the integral and inverting the kernel matrix. The subtraction point, s 0 is chosen to be near the Adler zero at leading order of χPT, s 0 4/3 m 2 π , and the subtraction constants g IL (s 0 ), which absorb the left hand cut contribution are the free parameters that are to be determined by fitting to the data.

III. NUMERICAL RESULTS
The η → 3π Dalitz plot distribution is conventionally expressed in terms of the variables x, y which are defined by For charge decay Q c = m η − 2 m + π − m 0 π and for the neutral decay Q n = m η − 3 m 0 π . Kinematics restrict the events to be contained within the unit disk x 2 + y 2 ≤ 1. The KLOE [18] and WASA-at-COSY [17] data were binned into 371 and 59 sectors of the unit disk, respectively (only bins that lie completely inside the physical I: Results of two-body (2b) and three-body (3b) fits to WASA-at-COSY [17] data, KLOE data [18] and the combined fit. For two-body fits we quote (g 2b IL (s0) ± ∆g 2b IL (s0))/g 2b 00 (s0), while when presenting results of threebody fit we quote (g 3b IL (s0)±∆g 3b IL (s0))/g 2b 00 (s0), where g 2b 00 (s0) is the central value obtained in the two-body fit with the same number of partial waves. We do the latter to illustrate the relative change in normalization between two-and three-body fits.
In the first step we fit the KLOE data alone. When only (I, L) = (0, 0), (1, 1) amplitudes are taken into account, we observe a significant reduction of χ 2 /d.o.f while moving from the "two-body" to the "three body" case. At the same time, when a complete set of S and P waves is incorporated, the χ 2 /d.o.f stabilizes at around 1.2-1.3 in both cases. In the second step, we combine the KLOE and WASA-at-COSY data. The results are in general very similar, showing the consistency of two different data sets. The results of the fit are shown in Fig. 1.

A. Matching to χPT and the Q-value
We note, that NLO χPT result depends on four low energy constants (LECs). These can be reduced to a single one L 3 = (−2.35 ± 0.37) × 10 −3 [50] if one employs Gell-Mann-Okubo constraint between meson masses and meson decay constants. This is not the case at NNLO where one has to deal with several unknown LEC's. Therefore, in our analysis we match the dispersive amplitudes with NLO χPT near the Adler zero. Note, that we match single variable partial wave amplitudes a IL (s) to χPT and not the full amplitude A C (s, t, u) along the lines s = t or t = u. This procedure should be equivalent, since a χP T IL (s) possess Adler zeros as well. In order to perform the matching at t = u we would need to make an additional analytic continuation of our results. In Fig. 2 we show our results of a combined fit to KLOE and WASAat-COSY with a fixed overall normalization to NLO χPT near the Adler zero in the region s = (0, 10 m 2 π ) along the lines s = t or t = u. The updated Q-value is which should be compared to the result of [28] Q = 21.4 ± 1.1 (the fit to WASA-at-COSY data only) and Q = 21.7 ± 1.1 (the fit to KLOE data only). Note, that the obtained Q-value is consistent with the latest (N f = 2 + 1 + 1) lattice computations Q = 22.2 ± 1.8 [51].
There are several challenges in the accurate determination of the Q-value. The first one comes from the elastic ππ scattering amplitudes, which are available from the Roy equation analysis [52] and implies the error ∆Q ππ = 0.25. Second uncertainty is due to experimental η → π + π − π 0 decay width, which serves as an input in our analysis. Its value increased by more than 3σ over the last thirty years, resulting in the current PDG value Γ η→π + π − π 0 = 296 ± 16 eV [1]. This error propagates to ∆Q Γ = 0.29. Third source of uncertainty is the experimental data on Dalitz plot itself, which thanks to the recent high-statistical analyses has improved significantly. Its contribution to the Q-value is < 10% of the size of error bars coming from the ππ amplitudes and therefore we included it in ∆Q ππ . Another uncertainty comes from matching to χPT amplitude. The error associated with L 3 LEC is very small and therefore the resulting error bar in our previous analysis was ∆Q total = 0.4 [28]. That error was dominated by the experimental error bars and therefore should be viewed as a lower bound of the full error. We note, however, that the Q-value determination is very sensitive to the matching to NLO amplitude. Though the region around the Adler zero is supposed to be stable against contributions from higher orders in the chiral expansion, we cannot completely exclude them. Assuming conservatively an additional error of 10% on NLO χPT amplitude, gives ∆Q match = 1.08 and the total ∆Q total = 1.1 quoted in Eq. (13).

IV. CONCLUSIONS
In this work we revisited our previous dispersive analysis [28] of the η → 3π decay in light of the new KLOE [18] data. Within our unitary model we established a unified description of charged and neutral decay modes. The method is based on Khuri-Treiman equation which is consistent with elastic unitarity, analyticity and crossing symmetry. Using the input from the ππ amplitude, the Khuri-Treiman equation was solved using Pasquier inversion technique. This allowed to establish a significant reduction of the unknown parameters compared to a more straightforward Omnès solution. However, the price is the treatment of the left-hand cuts, which is in general not known. We assume, that the unitarity in the physical region, where it can be constrained by the data, plays the key role and does not depend on an accurate form of the The amplitude along the lines s = u and s = t with a comparison to NLO ChPT. The relation between M (s, t, u) and A C (s, t, u) is given in [28]. The solid thick and thin vertical lines correspond to the physical region and the region where we calculated our amplitudes, respectively. As explained in the text, in order to compute the amplitudes beyond that region, a proper analytical continuation is required within Pasquier inversion technique.
unphysical left-hand cuts. The latter we absorbed in the subtraction constants [43]. With these model assumptions we were able to describe the data from KLOE [18] and WASA-at-COSY [17] with a minimal number fitting parameters.
Since the experimental data on η → 3π Dalitz plot is very precise now, the main experimental uncertainties come from I = 0 two-pion scattering amplitudes and the decay width Γ η→π + π − π 0 . Improving them are relevant for further Q-value and α determinations.
After submission of our manuscript an improved dispersive analysis based on Omnès functions was announced in [53]. The new Q-value is Q = 22.0 ± 0.7 which is consistent with our estimate.
The codes employed to compute the partial wave amplitudes and the Dalitz plot distribution are available for downloading as well as in an interactive form online at the Joint Physics Analysis Center (JPAC) webpage [54].