$\eta^{\prime}\to\eta\pi\pi$ decays in unitarized resonance chiral theory

We study the hadronic $\eta^{\prime}\to\eta\pi\pi$ decays within the framework of $U(3)_{L}\otimes U(3)_{R}$ Chiral Perturbation Theory including resonance states and the complete one-loop corrections. The amplitude is projected in partial waves and unitarized by means of the $N/D$ method resumming both the important $S$-and $D$-wave $\pi\pi$ and the subleading $S$-wave $\pi\eta$ final-state interactions. The participating scalar multiplet mass and coupling strengths are determined from fits to the Dalitz plot experimental data recently released by the A2 collaboration. As a byproduct of our analysis, the associated Dalitz-plot slope parameters are found to be $a=-0.072(7)_{\rm{stat}}(8)_{\rm{syst}}\,, b=-0.052(1)_{\rm{stat}}(2)_{\rm{syst}}\,, d=-0.051(8)_{\rm{stat}}(6)_{\rm{syst}}$, which lie in the ballpark of the current experimental and theoretical determinations.


Introduction
Different to Quantum Electrodynamics, a perturbative expansion in terms of the strong coupling cannot be applied to describe QCD processes at low-energies because the coupling strength becomes very large and therefore invalidates such an expansion. A well-known and particularly successful approach to overcome this limitation is Chiral Perturbation Theory (ChPT) [1], the low-energy effective field theory of QCD. ChPT is described in terms of eight pseudo-Goldstone bosons associated to the spontaneous chiral symmetry breaking SU (3) L ⊗ SU (3) R → SU (3) V exhibited by QCD i.e. three pions π ±,0 , four kaons K ± , K 0 andK 0 , and the η. The theory is constructed by performing a double perturbative expansion, in momenta, p 2 , and quark masses, m q , over Λ χ ∼ m ρ ∼ 1 GeV 1 . ChPT has been successfully applied for describing numerous processes involving pions and kaons but much less for the η. Actually, the η entering ChPT is not the physical one but rather a part of it corresponding to the octet 2 . In reality, the η meson has a second component, coming from the pseudoscalar singlet η 1 , which is not systematically included in ChPT due to the emergence of an anomaly. Indeed the U (1) A symmetry is broken (even in the massless case) by the quantum dynamics of QCD itself, preventing the η 1 to be the ninth Goldstone boson. This makes the η too heavy to be included as the ninth pseudo-Goldstone boson. However, in the limit of the number of colours becoming large, the "large-N C limit", the axial anomaly vanishes and the η 1 can be integrated to the Goldstone bosons [2,3,4]. In this limit, the (inverse) number of colors 1/N C is included in the power counting scheme as δ ≡ {(p/Λ χ ) 2 , m q /Λ χ , 1/N C } leading to a combined triple expansion in δ ∼ p 2 /Λ 2 χ ∼ m q /Λ χ ∼ 1/N C . Moreover, the SU (3) L ⊗ SU (3) R symmetry is enlarged to U (3) L ⊗ U (3) R and, the pseudoscalar octet and singlet states η 8 and η 1 mix allowing for a reasonable dynamical description of the physical η and η mesons.
While the convergence of the SU (3) ChPT perturbative expansion is restricted to lowenergies i.e. when the energy available in the process is below the mass of the first resonance (i.e. the mass of the ρ(770)), the η meson in U (3) is heavier (m η ∼ 958 MeV) than some resonances. To describe processes involving the η meson the ChPT framework has therefore to be enlarged to include explicitly such resonances. This is the avenue of Resonance Chiral Theory (RχT) [5]. In this theory, the interactions of the pseudoscalar mesons are governed by resonance exchanges. For these reasons, predicting observables that include η and η mesons is, typically, more difficult than for pions and kaons.
This article is structured as follows. In section 2, we define the kinematics of the process, introduce the Dalitz-plot parametrisation and discuss the current status of the associated parameters. The relevant Lagrangian is given in section 3. The structure of the decay amplitude is addressed in section 4 while section 5 is devoted to its unitarization. In section 6 we present the results of our fit to the experimental data for the π 0 π 0 mode from A2. Different fits are performed. They are organized according to their increasing fulfilment of unitarity. In each of these fits, the mass and the couplings of the scalar resonances are determined as well as the Dalitz plot parameters. We start our study considering ChPT including resonances and one-loop corrections. In a second step, we show the importance of the ππ S-and D-wave rescattering effects that nicely accommodate the π + π − cusp effect seen for the first time in η → ηπ 0 π 0 by the A2 collaboration. We then explore first the individual effect of the πη final-state interactions, anticipated to be small, before presenting our central description of the process including both ππ and πη rescattering effects. Our results are obtained from a fit to the A2 data. Our analysis enables us to extract some information about the I = 1 πη scattering phase shift within the allowed physical decay region. Moreover we predict the Dalitz-plot parameters and distribution of the π + π − decay channel that are found to be in excellent agreement with the BESIII experimental data. Finally, our conclusions are presented in section 7.

Kinematics and Dalitz-plot parametrisation
Let us consider the amplitude for the η (p η ) → η(p η )π(p 1 )π(p 2 ) decay, M(s, t, u). It is given in terms of the Mandelstam variables s = (p η − p η ) 2 = (p 1 + p 2 ) 2 , t = (p η − p 1 ) 2 = (p η + p 2 ) 2 , u = (p η − p 2 ) 2 = (p η + p 1 ) 2 , (2.1) which fulfill the relation The partial decay rate reads [41] Γ (η → ηππ) = 1 256π 3 m 3 η N ds dt |M(s, t, u)| 2 , (2.3) where N accounts for the number of identical particles in the final state; N = 1 for the charged and N = 2 for the neutral decay modes, respectively. The boundaries of the physical decay region in t lie within [t min (s), t max (s)] with t max/min (s) = 1 2   m 2 η + m 2 η + 2m 2 π − s ± λ 1/2 (s, m 2 η , m 2 η )λ 1/2 (s, m 2 π , m 2 π ) s   , (2.4) where λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz. The boundaries in s are given by However, the experimental measurements are often given as a power expansion in terms of the so-called Dalitz variables X and Y . These two variables are defined by where T π 1,2 and T η are the kinetic energies of the mesons in the η rest frame: and The decay width Eq. (2.3) can therefore also be written as: where now the integration boundaries are given by and The function h(s) is defined as (2.11) Fig. 1 shows the boundaries of the Dalitz plot in m 2 ππ ≡ s and m 2 πη ≡ t, u, the invariant masses (left panel) and in the Dalitz variables X and Y (right panel). At the ππ threshold, m 2 πη ∼ 0.59 GeV 2 , the two pions move together in the same direction with equal velocities and the η moves in opposite direction (red point). At m 2 ππ ∼ 0.11 GeV 2 , the range of m 2 πη increases going from ∼ 0.47 GeV 2 to ∼ 0.67 GeV 2 . In this last point, one pion is at rest and the other one moves in opposite direction of the η (blue diamond). At m 2 πη ∼ 0.54 GeV 2 , the allowed range of m 2 ππ values is large and reaches energies close to the region of influence of the σ meson. When m 2 ππ ∼ 0.16 GeV 2 , the η is at rest and the two pions move back-to-back (green triangle). We can anticipate that the region around this point will contain the largest number of events of the Dalitz plot decay distribution. Finally, at the πη threshold, m 2 ππ ∼ 0.12 GeV 2 , the η and one pion move in one direction with equal velocities and the other pion in the opposite direction (orange square).
The Dalitz plot parametrisation for η → ηππ decays is obtained by expanding the squared of the decay amplitude in powers of X and Y around the center of the Dalitz plot   Table 1 contains the current state-of-the-art Dalitz-plot parameters extracted from measurements together with their theoretical estimates. We can see large discrepancies 4 An alternative parameterization would be the so-called linear expansion where α is complex. A comparison with the parameterization of Eq. (2.12) gives a = 2Re(α) and b = Re(α) 2 + Im(α) 2 . The two parameterizations are equivalent if b > a 2 /4.
# events GAMS-4π [43] −0.067 (16) [42] −0.127 (16) [39] is in agreement with VES result. This is not surprising since this data set was fitted in their analysis. In Ref. [35] the parameter a is not predicted but rather fixed from an average of the VES and GAMS-4π values. From fits to the 2011 BESIII and VES data, the recent dispersive analysis of Ref. [40] obtains values ranging from −0.041 to −0.148, which agree with the corresponding measured values reported by these two collaborations. The parameter b determined by VES is ∼ 1σ away from the measurements of GAMS-4π, BESIII and A2. A deviation of ∼ (1-2)σ is also seen with respect to the analyses of Refs. [35,39]. However, this discrepancy disappears in the analysis of Ref. [40]. Both recent experimental measurements and theoretical predictions seem to indicate an unambiguously negative value for b in clear disagreement with a vanishing b obtained in Refs. [37,47]. Regarding the parameter c, the symmetry of the wave function forbids such a term in the neutral channel. In the charged channel the odd terms in X are forbidden since C-parity is conserved by strong interactions. So its value is predicted theoretically to be zero in the Standard Model. And the measured values of this parameter are consistent with zero. Finally, for the value of the parameter d, experimental results seem to favour the predictions of Refs. [35,40] with respect to those of Refs. [39,47]. In conclusion, from Table 1, we observe an inconsistent picture so far. The theoretical determinations of the Dalitz parameters are intricately linked to the data used to constrain the corresponding theories or models.

Formalism
The Large-N C ChPT is the effective field theory of QCD in the chiral and large-N C (number of colour) limits. Within this framework the singlet field η 1 , absent in SU (3) ChPT, becomes a new degree of freedom of the effective theory i.e. the ninth Goldstone boson associated with the spontaneous breaking of The theory is usually called U (3) ChPT and the same counting is assigned to the squared momenta p 2 , to the light quark masses m q and to the inverse of the number of colors 1/N C , giving rise to a combined triple expansion in δ ∼ p 2 ∼ m q ∼ 1/N C . With the use of the single power-counting parameter δ the expansion of the effective Lagrangian is given by [2] In this notation, the contributions L δ i are of order O(δ i ). For example, the leading and next-to-leading order Lagrangians, L δ 0 and L δ 1 , are of order O(δ 0 ) = O(1) and O(δ), respectively. One interesting feature of the combined power counting is that meson loop diagrams with vertices from L δ 0 count as O(δ 2 ). They are 1/N C suppressed with respect to the NLO diagrams from L δ 1 . Therefore at O(δ) only tree level diagrams from L δ 0 and L δ 1 need to be taken into account. The loop diagrams are higher order in the counting. In our analysis, we include the one-loop corrections for the first time for describing η → ηππ and work at O(δ 2 ) in order to match the high level of precision of the experimental measurements. We will study the impact of such inclusion. Contrary to SU (3) or U (3) ChPT, the inclusion of resonances into the description of the effective field theory spoils the power counting. This is why, while RChT includes ChPT at O(p 2 ) = O(δ 0 ), it does not include the next order either in the chiral or in the combined triple expansion scheme. It rather substitutes it by a Lagrangian accounting for the interactions between pseudoscalar mesons and resonances. In fact, resonance exchanges saturate the higher order local contributions of the effective expansion. As a consequence, the systematic effective field theory with a rigorous power counting scheme is lost. It is replaced by a model based on the large-N C limit as a guiding principle. This model allows in general for reasonable descriptions of QCD processes at low-energies.
The relevant Lagrangian for our work including scalar resonances is written as [5,48] where the dots denote operators with three or more resonance fields which we neglect in this analysis. The first term on the right-hand side of Eq. (3.2) is the chiral Lagrangian at leading order in U (3) ChPT, it is O(δ 0 ) and reads where the last term accounts for the U (1) A anomaly contribution to the pseudoscalar singlet η 1 of mass m 2 0 . The chiral building blocks are defined by F is the axial decay constant of the pseudo-Goldstone bosons in the chiral and Large-N C limits. s, p, v µ , a µ stand for external fields and the parameter B is related to the quark condensate 0|q i q j |0 = −F 2 Bδ ij . In the absence of external fields, i.e. v µ = a µ = p = 0, we have s = diag(m u , m d , m s ), with m q the light quark masses, encoding the explicit chiral symmetry breaking. The pseudo-Goldstone bosons are collected in the matrix 3) in terms of Φ, one obtains the kinetic term plus a tower of derivative interactions increasing in an even number of pseudoscalar fields. This together with the quadratic mass term plus additional interactions proportional to the quark masses gives (3.6) The mass term M in the previous expression induces a η 8 -η 1 mixing. At leading order, the physical η and η mass eigenstates are then obtained after diagonalising the mass matrix with the following orthogonal transformation: where the mixing angle θ stands for the degree of admixture. When higher order terms are to be considered, the transformations accounting for the mixing are more involved since not only the mass terms do mix but also the kinetic ones. In this work, we treat the η-η mixing as described in section III of Ref. [48] and refer the reader to this reference for more details and the expression of the higher order terms. In particular, we use the next-to-leading order expression given in Eq. (15) of Ref. [48] based on the one-angle scheme approximation 5 .
The second term in Eq. (3.2) corresponds to the interaction terms of two pseudo-Goldstone-bosons with one resonance and is given by The resonance state building blocks are (3.9) Similarly to Eq. (3.6), expanding the Lagrangian in Eq. (3.8) in terms of Φ we get where we have used χ = 2B 0 M and neglected other external fields (v µ = a µ = p = 0). Note that the interaction terms proportional to c d andc d enter only with derivatives while c m andc m are proportional to the quark masses. The third term in Eq. (3.2) contains the kinetic term, Finally, the last term in Eq. (3.2) is a local operator of O(δ), influencing only the singlet sector. It cannot be generated from the exchange of the scalar resonance discussed above. It is obtained by integrating out pseudoscalar resonances instead. It only involves pseudo-Goldstone bosons and reads where with a µ = (r µ − l µ )/2.

Structure of the decay amplitude
The calculation of the η → ηππ decay amplitude includes the diagrams depicted in Fig. 2 and can be gathered as M (2) is the amplitude at lowest order in ChPT. M Res represents the amplitude involving the exchange of resonances, M Loop the loop contributions, and M Λ the Λ term 6 . The lowest order ChPT contribution is constant where isospin symmetry m u = m d =m ≡ mu+m d 2 has been assumed. The origin of this term stems entirely from the interactions proportional to quark masses appearing in the last term of Eq. (3.6) i.e. no derivative interactions contribute at this order. It is thus chirally suppressed explaining the smallness of this contribution.
The (zero-width) resonance exchange contribution M Res reads Unitarity loop corrections to the decay amplitude can occur either in the s or in the t and u channels. In the s-channel, the meson pairs ππ, KK, ηη, ηη and η η enter within the loop while KK, πη and πη pairs contribute in the t and u-channels. Tadpole contributions with π and K also appear. The complete expression for the unitarity loop corrections M Loop , including the tadpoles can be found in appendix A. Finally, the term from L Λ is constant Moreover, there is a contribution from the mixing [48] with sin θ δ 0 [48] and . Although the mixing contribution is omitted all along the formulae presented in the next sections, we have included it in our analysis. Note that both the Λ and mixing terms, Eqs. (4.4) and (4.5), are constant and their contributions are found to be small. Indeed they are both proportional to m 2 π and hence chirally suppressed.

Partial waves and unitarisation of the amplitude
Within the ChPT framework described in sections 3 and 4, the loop contributions are unitary order-by-order in the perturbative expansion. Nearby the resonance region, however, the perturbative chiral amplitudes violate unitarity. This inherent limitation of the theory is addressed by using a unitarization procedure. For this, we rely on the N/D method. A detailed account of this method can be found in Refs. [48,55,61,62]. In the following, we recall the main features of this approach that are relevant for our analysis. Another method relying on the Khuri-Treiman framework has been followed in Ref. [40]. We start by writing the most general unitarity relation for η (p η ) → η(p η )π(p 1 )π(p 2 ) decay M η →n denotes the η → n decay amplitude and T n→ηππ the transition n → ηππ. If for simplicity we restrict n ≤ 3 then the scattering matrix element T n→ηππ contains pure 3body → ηππ contributions as well as two-body final-state interactions. The three-body final-state interactions are suppressed by power counting and phase-space compared to the two-body final-state interactions. In our analysis, we include only the dominant two-body final-state interactions. In a three-body decay, two-body final-state interactions can occur either by means of a rescattering where two out of the three final state particles rescatter an arbitrary number of times in each of the two-particle channels considering the third particle as a spectator or by interactions among one of the two rescattering particles together with the third spectating-particle. While the former will be fully accounted for in our study only portions of the later will be incorporated. In the following, we limit the sum over n in Eq. (5.1) to ππ and πη intermediate states. The unitarity condition for the η → ηππ decay in terms of the ππ and πη scattering amplitudes can be written as The symmetry factor is N = 2 in case of identical (ππ) and N = 1 for distinguishable (πη) particles. θ s,t,u stands for the center-of-mass scattering angle between the initial and intermediate state and θ s,t,u denotes the center-of-mass scattering angle between the intermediate and final state.
The decay and scattering amplitudes, M I and T I , can be decomposed in partial waves with definite isospin I and angular momentum J through where P J is the Legendre polynomial of the J th degree. Isospin conservation constrains the total isospin of the final state ππ and πη pairs to be I = 0 and I = 1, respectively. We limit our study to the S (J = 0) and D (J = 2) waves for the ππ scattering and to the S-wave for the πη scattering 7 . Inserting Eqs. (5.3) and (5.4) in Eq. (5.2) and integrating over the momentum using the relation the following unitarity relations for each partial-wave of the decay amplitude of definite isospin can be derived:

ππ final-state interactions
Let us consider first the unitarity relation in the s-channel for illustrating the N/D method applied to the η → ηππ decay. For a well-defined isospin I and angular momentum J from Eq. (5.6) we have A possible way to fulfill the previous equation is as follows. We write the partial wave associated to the perturbative decay amplitude, Eq. (4.1), as , (5.9) where "pert" stands for perturbative. The amplitudes with superscript (2) correspond to the tree-level amplitudes, Res to the resonance exchanges in the s, t and u channels, Loop to the loop contributions in the t and u channels and finally Λ represents the contribution from the Λ term. The function g ππ (s) entering the second term of Eq. (5.9) accounts for the discontinuity along the right-hand cut due to the two-pion intermediate states a ππ (µ) ≡ a ππ is a subtraction constant that is not directly determined by the unitarization procedure but should be fixed from elsewhere. g ππ (s) satisfies Img ππ (s) = −σ π (s)/16π [48].
It is related to the standard one-loop function given in Eq. (A.4) through where a ππ (µ) is an arbitrary subtraction constant. The basic idea of the N/D unitarisation method consists in collecting the left-and right-hand cuts in two different functions. A N/D representation of the decay amplitude, Eq. (5.9), can be obtained rewriting it as [54] contain, respectively, the corresponding perturbative calculation of the partial wave ππ scattering and the η → ηππ decay amplitudes. As explained previously, in t IJ ππ (s) , Eq. (5.13), the part with superscript (2) corresponds to the tree level amplitude, the one with superscript Res to the exchange of resonances in the s, t and u channels and the Loop one denotes the loop contributions in the t and u channels as well as the inelastic loop contributions in the s-channel. The corresponding diagrams are depicted in Fig. 3 recovering Eq. (5.9) up to higher orders. Note that a chiral expansion of Eq. (5.12) leads to + · · · . (5.14) π π π π π π π π S 8 , S 1 s-channel π π π π K , η , η , η We would like to emphasize that the function R(s) entering Eq. (5.13) does contain the left-hand cut (LHC) of the decay amplitude, perturbatively treated, but does not contain the ππ right-hand cut which is treated non-perturbatively using the function g ππ (s).
Thus, the absorptive part of the unitarised partial wave decay amplitude m IJ (s) as given in Eq. (5.12) satisfies s-channel unitarity 8 : 8 In principle, the unitarity relation given in Eq. (5.15) is valid up to the first inelastic threshold i.e. the KK threshold. However, there is a spurious contribution coming from the imaginary part of the t-and u-channel πη left-hand cut loops. They sit on the elastic region and induces a violation of unitarity. We have checked that this unitarity violation is numerically tiny and therefore acceptable for our purposes. There is, moreover, another source of unitarity violation coming from the unitarization method itself. Spurious singularities, as for example the ones given by the diagram below, η η π π π π π π are generated from the on-shell approximation within the N/D method by the t-and u-channel ππ loops entering the ππ scattering amplitude. In fact, this a drawback of the unitarization method (see Refs. [49,56] for more details about this pathology). This violation of unitarity turns out to be acceptable since these effects are usually found to be numerically small. See section 6.2 for a discussion of the size of these constributions in the present case. .
The unitarized η → ηππ I = 0 decay amplitude written in terms of the S-and D-waves is with t max/min defined in Eq. (2.4) and where θ s is the angle of p π with respect to p η in the rest frame of the pion pair. It is given by .

πη final-state interactions
By analogy with Eq. (5.12), the N/D representation of the decay amplitude can be written as where the πη final-state interactions in the t-and u-channels have been resummed. We have contains the perturbative calculation of the partial wave πη scattering and m IJ η →ηππ (s) (2)+Res+Loop+Λ the η → ηππ decay amplitude one. The diagrammatic structure of the πη → πη scattering amplitude resembles the ones shown in Fig. 3 for ππ. The notation here follows the convention adopted for ππ see Eq. (5.13) with ππ replaced by the πη-system.
The two-particle discontinuity along the right-hand cut due to the πη intermediate states reads and analogously for g πη (u) with t ↔ u. However, we slightly modify the expression of Eq. (5.20) and construct the amplitude such that the perturbative terms of the decay amplitude are kept. They are supplemented by the inclusion of the S-wave πη final-state interactions. In this way, we are able to quantify the importance of the πη rescattering contribution with respect to the perturbative calculation. Since Eq. (5.20) generates the S-wave projection of the lowest-order, the resonance exchanges, the loop contributions and the Λ term 9 , we need to remove these contributions and add by hand the term M(s, t, u) (2)+Res+Loop+Λ , see section 4. Thus, the unitarized decay amplitude for η → ηππ taking into account the I = 1 S-wave πη final-state interactions can be written as Note that M(s, t, u) (2)+Res+Loop+Λ is not projected into partial waves and the last two terms are introduced to avoid double counting. The m 10 η →ηππ partial waves are derived through where , (5.26) with ∆ P Q = m 2 P − m 2 Q . In Eq. (5.25), θ t is the angle of p η with respect to p π in the πη center-of-mass frame given by .
Analogous expressions are valid in the u-channel with the replacements t ↔ u and cos θ t ↔ − cos θ u . The required πη-scattering amplitude and the calculation of the corresponding S-wave projection t 10 πη (s) (2)+Res+Loop+Λ are given in appendix C.

Fits to experimental data
We relate the theoretical expression for the differential decay rate of η → ηπ 0 π 0 , Eq. (2.8), to the Dalitz distribution of the measured number of events through For our study, we analyze the acceptance corrected η → ηπ 0 π 0 Dalitz distribution recently released by the A2 collaboration [45]. The factor of 2 accounts for the fact that the data is given for half of the (symmetric) Dalitz distribution. N events is the total number of events for the considered process. Γ η is the total decay width of the η meson. ∆X and ∆Y are the bin width of the X and Y variables.B(η → ηπ 0 π 0 ) ≡B is a normalisation constant that, for a perfect description of the spectrum, would be equal to the corresponding branching fraction. For our analysis, we fix this normalisation to the PDG reported valueB = 22.8(8)% [41] 10 . Two analyses of the same data set called analysis I and analysis II have been performed in Ref. [45]. They correspond to different analysis frameworks and to different selections of data samples. The corresponding efficiency corrected numbers of events (N events ) are 463066 for analysis I and 473044 for analysis II. The width of the bins is ∆X = ∆Y = 0.10 MeV. The central value of the results of our analysis are obtained by considering the data set of analysis I. The results obtained with the data set of analysis II is used to assess the systematic uncertainties of our fit results presented in section 6.3. The χ 2 function we minimize is 2) where N exp XY is the experimental number of events and σ exp XY the corresponding uncertainties in the XY -th bin. The number of data points to be fitted is 200.
Our fitting strategy is organized in a bottom-up approach guided by step-by-step implementation of two-body unitarity. The free parameters to fit are M S 8 , M S 1 , M a 0 , c m ,c m , c d andc d . However, in order to reduce the number of free parameters we invoke the large-N C relations for the couplings and masses of the octet and singlet. We [5]. For the η-η mixing angle we take θ = −13.3(5) • [63] 11 while we use Λ 2 = −0.37 [54] and set the regularization scale to µ = 700 MeV. Note that by imposing the large-N C constraints we are introducing large correlations between the fit parameters. This implies that different values for the couplings and masses can lead to a similar fit quality underestimating the true statistical errors. Therefore the fit results presented below should be taken with a word of caution.

ChPT fits including resonances and one-loop corrections
We start by fitting Eq. (6.1) to A2 data [45] with the amplitude described in section 4. For our first fit we impose the restriction c d = c m for the couplings which comes from the requirement of the Kπ scalar form factor to vanish at high energies [64]. The resulting fit parameters take the values The fitted parameters are particularly strongly correlated in this case leading to a large error for the coupling c m . This is not a surprise since c m always enters with m 2 π in the amplitude, and hence is chirally suppressed, indicating that its influence is small. This is in agreement with previous estimates of this coupling suffering from a large uncertainty. For example, c m = 31.5 +19.5 −22.5 MeV in [48], c m = 15(30) MeV in [61] and c m = 80(21) MeV in [52]. In order to alleviate this correlation, we also consider fits where we fix the couplings c m andc m . For instance, we first fix them to c m = 41.1(1) MeV andc m = 18.9(9) MeV [56] using results obtained from meson-meson scattering (see Refs. [48,54]  We have also tried fits having all couplings i.e. c m ,c m , c d andc d as free parameters but these fits are unstable because there are too many free parameters to fit. In Fig. 4, we provide a graphical account of the ratio of the differential decay width distributions of η → ηπ 0 π 0 as a function of X, Y, m π 0 π 0 and m π 0 η over the phase-space obtained from the fit results of Eq. (6.3). In order to compare with the experimental data, both expressions are normalized such that the individual integrated branching ratio is 1. The corresponding normalized amplitudes are denoted asM andφ in the figure. A cusp effect at the π + π − mass threshold is neatly visible in the data (see top-right and bottomleft panels of the figure). This is the first cusp structure observed in η → ηπ 0 π 0 and is not accommodated by the theoretical description given here. In the next section 6.2 we will improve the theoretical description in order to obtain a better agreement with the experimental data and to try to describe the cusp effect.
In Fig. 5 we display graphically the different contributions entering the decay amplitude as a function of the m π 0 π 0 invariant mass distribution. The hierarchy between the resonance exchange and loop contributions is shown. From the top-left panel, we observe that the lowest-order contribution is tiny while the decay is largely dominated by the resonance exchanges with the loop contributions interfering destructively. The integrated branching ratio associated to these curves is of 0.6% for the lowest-order (blue dotted curve), 28.2% for the lowest-order plus resonance exchanges (red dashed curve) and of 22.5% for the lowest-order including resonance exchanges and loop contributions (black solid curve).   Finally, we would also like to provide an estimate of the chiral couplings L 5 and L 8 as well as the sum 3L 2 + L 3 as output of our fits. Assuming resonance saturation for the order p 4 ChPT couplings constants [5], the following relations can be derived  With the results of the fit Eq. (6.7), we get: 3L 2 + L 3 = 0.51 · 10 −3 , L 5 = 1.03 · 10 −3 , L 8 = 0.51 · 10 −3 , (6.15) Considering the results of the fit Eq. (6.11), we obtain The values of Eqs. (6.14) and (6.15) are in reasonable agreement with the chiral couplings determinations at O(p 4 ) (see Ref. [51] for a recent review) while the values for L 5,8 in Eq. (6.16) disagree. However, the first relation in Eq. (6.13) is not well fulfilled for most of the values of L 2 and L 3 obtained at O(p 6 ). For example, the most recent analysis of K 4 decays gives L r 2 = 0.63(13) · 10 −3 and L r 3 = −2.63(46) · 10 −3 [50] leading to a negative value for the sum 3L 2 + L 3 . So the left-hand side of the first equality of Eq. (6.16) is negative while the right-hand side of the equation is definite positive. This inconsistency shows that resonance saturation of low-energy constants by scalar resonances should be taken with a word of caution.

Fits including individual ππ and πη final-state interactions
In order to improve our fits to the data, we unitarize the parameterization of the decay amplitude. As derived in Eq.  (14) , (6.18) for the associated Dalitz-plot parameters. Note that, with respect to the previous section, there is one more free parameter to fit since we have the subtraction constant, Eq. (5.11) to determine. Due to the important correlations between the fit parameters, the statistical uncertainties have increased while the total χ 2 /dof is slightly improved. Comparing these results to the ones of Eq. (6.4) obtained within ChPT including resonance exchanges and one-loop corrections without considering any final-state interactions, shows that the inclusion of the ππ final-state interactions has a large effects on the Dalitz plot parameters.
In particular, the a and b parameters associated to powers of Y , have been substantially shifted downwards while the change on the d parameter related to X 2 is slightly less severe. We would like to note that while the S-wave only affects the parameters a and b associated to powers of Y , the D-wave affects the variable X 2 but also the determination of a and b.
For example, if we consider only the S-wave for the ππ final-state interaction, we obtain a = −0.094 (7) , b = −0.034(2) , (6.19) instead of the results of Eq. (6.18). Comparing to Eq. (6.4) a seems unaffected by the inclusion of the S-wave ππ final-state interaction while b is clearly moved down. We therefore conclude that while a precise determination of the b parameter requires to take into account the the S-wave ππ final-state interaction the determination of the a and d parameters are dominated by the D-wave.
We shall now return to the discussion on the spurious singularities introducing unitarity violations, see section 5.1. According to the Watson's theorem [65], the phase of the η → ηπ 0 π 0 decay amplitude equals the ππ scattering phase-shift in the elastic region. On figure 6 we compare the S-wave phase of the decay amplitude (blue line) to the ππ scattering phase-shift (red dashed line) using the results of the fits Eq. (6.17). The difference between the two phases in the elastic region is due to such spurious contributions. This implies that unitarity is no longer fulfilled exactly. However, this effect is numerically quite small (for the D-wave this effect is completely negligible) and acceptable for our analysis. We have checked that by removing these spurious contributions (i.e. removing the t-and u-channel from the ππ loop contributions) entering ππ scattering, unitarity is restored. The changes in the fits are negligible.   In this case we find for the chiral couplings: for the fit Eq. (6.11).  table 2 to the ones given in section 6.1 obtained using ChPT including resonances and one-loop corrections without resumming the final-state interaction effects, we observe that the inclusion of the πη rescattering effects has small effects on the result of the fits. We therefore conclude that πη rescattering effects are small as previously observed, see Ref. [58]. Note that Fit B allowing the coupling c m to float carries a large error bar for the same reasons as discussed in the previous section.
As in section 6.1, we compare the fit results to the experimental data in Fig. 7. The black solid curve corresponds to the fit results Eq. (6.17) where the ππ final-state interactions have been taken into account. The gray dashed curve represents the resulting amplitude obtained from Fit A of table 2 for which the πη final-state interactions have been resummed. Contrary to the fit results shown in Fig. 4, the cusp structure at the π + π − mass threshold (see top-right and bottom-left panels of the figure) is now nicely accounted for within our description. This is possible after the inclusion of the ππ unitarization 13 . Note that the theory describes much better the data once ππ rescattering effects have been included compared to πη 14 . This is also in part reflected in the corresponding χ 2 dof .

Fits including both ππ and πη final-state interactions
This section contains our central results for the corresponding fit parameters as well as for the associated Dalitz-plot slope parameters. They are obtained by building a representation for the amplitude that takes into account the ππ and πη final-state interactions simultaneously as described in section 5.3. In the previous sections, we presented fits to the analysis I of the A2 data set. Here we fit our theoretical representation to analysis I and II of the A2 data set allowing us to include a systematic uncertainty on our results. The systematic uncertainty is taken to be the difference between the results of the fits to the two different experimental analyses. The fit results are collected in table 3 using two different settings and with fixing the πη subtraction constant to a πη = 2.0 +3.1 −3.4 [48]. i) In Fit 1 we impose two conditionsc d,m = c d,m / √ 3 and c d = c m ; ii) In Fit 2 we fix c m = 41.1(1) MeV andc m = 18.9(9) MeV [56] and imposec d = c d / √ 3. Contrary to what was done in the previous sections, we report here not only the results for the Dalitz plot parameters a, b and d but also the results for the higher order ones, κ 03 , κ 21 and κ 22 . They are found to be very small as reported in previous theoretical analyses [35,40].
Note that the values found for the Dalitz plot parameters do not change much with the different fit scenarios and are very similar to the ones obtained in section 6.2 where only the ππ rescattering effects were taken into account. This is expected since we saw that the ππ rescattering dominates the final-state interactions. The stability of our fit results makes us very confident in the robustness of the results.
In order to illustrate the overall effects of the D-wave ππ final-state interactions, we have also performed fits to the Dalitz plot experimental distribution without the D-wave contribution. The resulting fit results are gathered in   Table 3: Results for the parameters of the fits together with their associated Dalitz parameters for two different fit scenarios and two different analyses (analysis I and II) of the A2 data set. Masses and coupling are given in MeV. The first error is the statistical uncertainty coming from the statistical uncertainties on the data, the second error is the systematic uncertainty coming from the uncertainty on the subtraction constant a πη . See main text for details.       In the following, we study the dependence of the Dalitz parameters with respect to the numerical values of the mass and couplings of the participating scalar multiplets. For this exercise, we take different values for the mass and couplings from the literature and make some "crude" predictions. The resulting estimates are gathered in table 5 where we have used the constraintsc d = c d / √ 3 andc m = c m / √ 3, and fixed a ππ = 0.76 from table 3. These results show that the Dalitz plot parameters are sensitive mostly to the values of c d and M S . The variation of these parameters has also a sizeable impact on the predicted branching ratio. Out of the five predictions shown in this table, the results given in the last column are the most realistic ones.
In Fig. 8 we compare the experimental data to the results corresponding to Fit 1 of Parameter Predictions Constraint Ref. [5] Ref. [53] Ref. [53] Ref. [  We observe that the representation of the amplitude obtained from the fit results of analysis I (black solid curve) practically overlaps with the one coming from the fit results of analysis II (gray dotted curve). The representation of the amplitude built from the results of our fits successfully describes the experimental data including the cusp effect at the π + π − threshold. Moreover, in order to compare the decay amplitude with the Dalitz plot experimental measurements, we compute its square in terms of the Dalitz variables X and Y inside the physical decay region. The shape of the Dalitz distribution, normalized to 1 in the center of the Dalitz plot, is displayed in Fig. 9 for the ChPT results presented in Eq.  [45,46]. It shows that the Dalitz distribution is more populated when the pions go back-to-back (cf. Fig. 1). In order to further illustrate the strong effects of all finalstate interactions on the Dalitz plot distribution, on the bottom left panel of Fig. 9 we plot the quantity |M (X, Y ) Full | 2 divided by the same quantity before the unitarization, corresponding to |M (X, Y ) ChPT+Res+Loop | 2 . Clearly, the effects of the unitarization of the amplitude, dominated by the ππ rescattering, are very important in the upper central region of the distribution. Finally, we study the region of the Dalitz plot influenced by the effects of the D-wave ππ final-state interactions. The answer is given on the bottom right panel of Fig. 9 where we show the quantity |M (X, Y ) Full | 2 divided by |M (X, Y ) D−wave=0 | 2 . |M (X, Y ) D−wave=0 | 2 corresponds to |M (X, Y ) Full | 2 with the D-wave ππ final-state interactions effects set to zero. We can see that the D-wave effects also appear on the upper central region of the Dalitz plot.       table 3 (black solid and gray dotted curves for the analyses I and II data sets, respectively). They are obtained after resumming both the ππ and the πη final-state interactions. Data is taken from Ref. [45].
The central values of our final results for the Dalitz-plot parameters associated to the η → ηπ 0 π 0 decay correspond to Fit 1 of table 3 for the analysis I of A2 data. To assess the "experimental" systematic uncertainty, we take the largest variation of the central values with respect to the results considering the analysis II data set of the same table. We add this uncertainty to the systematic uncertainty coming from the subtraction constant a πη in quadrature and we obtain Although our dedicated analysis shows that the πη rescattering effects are small we can still extract some information about the I = 1 πη phase shift as a byproduct of our study. In Fig. 10  Before concluding, in Fig. 11, we compare our results on the π 0 π 0 mode to the BESIII measurement [46]. Contrary to A2, the BESIII experimental data is not yet publicly available, so we have extracted the data points from figure 7 of Ref. [46] for the comparison. Our prediction is displayed in Fig. 11. It is in very good agreement with the measured data. To show this we computed the χ 2 /dof using the results of Fit 1 of table 3 and the BESIII data. We obtain χ 2 /dof= 101.5/95 ∼ 1.07. Note that contrary to A2 no statistically significant evidence for a cusp at the π + π − threshold is observed.
Using our representation of the amplitude using the fits to the data on the π 0 π 0 mode from the A2 collaboration, we can predict the Decay rate distribution in the charged channel (π + π − ). To predict the Dalitz-plot parameters of the π + π − decay mode, one should This work: our prediction Figure 11: Differential decay rate distribution for η → ηπ 0 π 0 divided by the phase-space, both individually normalized, associated to the resulting parameters of Fit 1 of table 3 as compared with the BESIII experimental data [46]. consider all possible sources of isospin breaking. In our framework, isospin breaking effects mostly affect the Dalitz variables X and Y if the charged pion mass is used in Eqs. (2.9) and (2.10). In Ref. [40] relations between the Dalitz parameters in the charged and the neutral decay modes have been derived: where the superscripts c and n denote the associated parameters in the charged and neutral systems, respectively, and with ε iso ∼ 4.7% [40]. Following this prescription, our estimates for the Dalitz parameters in the charged channel reads a = −0.065(7) stat (8) syst , b = −0.048(1) stat (2) syst , d = −0.045(7) stat (5) syst . (6.31) Comparing the above results with the most recent experimental determination of these parameters in the charged system released by BESIII in 2017 [46], a = −0.056(4) stat (3) syst , b = −0.049(6) stat (6) syst , d = −0.063(4) stat (4) syst , we observe that our prediction for b is in excellent agreement while a and d are found to be 1σ and 2σ away, respectively.
Finally, our results given in Eqs. (6.27) and (6.31) for the neutral and charged decays modes, respectively, are graphically compared to previous experimental and theoretical determinations in Fig. 12 Figure 12: Comparison of experimental (•) and theoretical ( ) determinations of the associated Dalitz-plot slope parameters for η → ηππ (cf. table 1). Our results ( ) correspond to Eqs. (6.27) and (6.31) for the π 0 π 0 and π + π − modes, respectively. Only the statistical uncertainty is shown.

Conclusions
Recent measurements of the η-η system have reached unprecedented precision placing new demands on the accuracy of the corresponding theoretical description. The η → ηππ decays represent a good laboratory to test any extension of SU (3) Chiral Perturbation Theory, the effective field theory of QCD, which has proven to be very successful in describing pion and kaon processes. In this work, we have analyzed the η → ηππ transition within U (3) ChPT at one-loop including scalar resonance states as degrees of freedom.
The corresponding amplitude has been unitarized using the N/D method. Our treatment accounts for simultaneous ππ and πη final-state interaction effects. This parametrization has been fitted to the recently released A2 collaboration data on the η → ηπ 0 π 0 channel and very good agreement has been achieved. The results of the fit show that the Dalitz plot parameter b is shifted downwards compared to U (3) ChPT predictions. We demonstrate that this is attributed to the S-wave resummation of the ππ final-state interactions. Moreover, to match the A2 experimental accuracy requires the inclusion of the D-wave contribution. On the contrary, the S-wave πη rescattering has shown to be small in agreement with previous studies.
To further improve the description of the rescattering effects, one can consider a more sophisticated unitarization procedure in coupled channels including inelastic KK scattering. We postpone it for a future analysis when new measurements, e.g. by GlueX experiment, become available.
In summary, from our analysis we extract the following Dalitz-plot parameters a = −0.072(7) stat (8) syst , b = −0.052(1) stat (2) syst , d = −0.051(8) stat (6) syst . Using these results, we are able to make predictions for the charged channel. These predictions are found to be in very good agreement with the BES-III measurements of this channel. Moreover, we were able to extract some information on the I = 1 πη phase shift at low energy.
The theoretical framework developed here should be suitable for precision analyses of future experimental data.