Patterns of $C$- and $CP$-violation in hadronic $\eta$ and $\eta'$ three-body decays

We construct hadronic amplitudes for the three-body decays $\eta^{(\prime)}\to\pi^+\pi^-\pi^0$ and $\eta'\to\eta\pi^+\pi^-$ in a non-perturbative fashion, allowing for $C$- and $CP$-violating asymmetries in the $\pi^+\pi^-$ distributions. These amplitudes are consistent with the constraints of analyticity and unitarity. We find that the currently most accurate Dalitz-plot distributions taken by the KLOE-2 and BESIII collaborations confine the patterns of these asymmetries to a relative per mille level. Our dispersive representation allows us to extract the individual coupling strengths of the $C$- and $CP$-violating contributions arising from effective isoscalar and isotensor operators in $\eta^{(\prime)}\to\pi^+\pi^-\pi^0$ and an effective isovector operator in $\eta'\to\eta\pi^+\pi^-$, while the strongly different sensitivities to these operators can be understood from chiral power counting arguments.


Introduction
The idea that the conservation of discrete symmetries in the strong interactions does not have to be manifest first rose to prominence in the year 1950 with the work of Purcell and Ramsey [1], who proposed the violation of P and CP in the decay η → 2π. Later, this idea was theoretically realized in a P -and CP -odd operator of dimension four in QCD, which is well known as the θ-term. The latter induces, amongst others, an electric dipole moment (EDM) of the neutron. Rigorous experimental limits on EDMs imply corresponding theoretical limits on η → 2π [2][3][4], a link that can even be established without recourse to the θ-term as the fundamental mechanism [5][6][7][8][9][10]. Accordingly, no measurement so far could find evidence for this process, which is probably beyond experimental reach for the foreseeable future. Given the dearth of experimental evidence for sources of CP -violation beyond the Cabibbo-Kobayashi-Maskawa mechanism in the weak interactions of the Standard Model (SM), it is worthwhile to investigate another category of CP -violating operators that has gained much less attention so far: T -odd and P -even (TOPE) interactions, which in addition violate C according to the CP T theorem. C-violating effective operators have been discussed in the literature to some extent [11][12][13][14][15], but explicit links to hadronic processes have barely been established. In the Standard Model effective field theory [16,17], they only contribute starting at dimension 8 [18,19]. 1 As TOPE forces cannot be mediated by π 0 exchange [20], they can only contribute via short-range nuclear forces, and are therefore far less constrained from nuclear physics. Suitable candidates to investigate these kinds of operators are certain decays of the η (′) mesons, which are eigenstates of C. These allow us to investigate TOPE forces in the absence of the weak interaction, such that the observation of a corresponding C-violating η (′) decay would automatically indicate physics beyond the Standard Model (BSM).
Studying the charge asymmetry of the η → π + π − π 0 Dalitz-plot distribution offers an ideal stage in the search for such BSM physics. As pointed out in Ref. [21], in contrast to other C-violating processes such as η (′) → 3γ, η (′) → π 0 γ * , etc., the breaking of mirror symmetry in η → π + π − π 0 is linear in these BSM operators, as it is generated through interference with the SM mechanism. For an overview of C-and CP -violating processes in the η and η ′ sector we suggest Ref. [10]. The simplest observable that can be probed experimentally is the left-right asymmetry A LR that compares the two halves of the Dalitz-plot distribution divided along the π + ↔ π − line of reflection [22]. It is also possible to construct more sophisticated quadrant and sextant asymmetry parameters A Q and A S that allow us to disentangle the contributions of the BSM ∆I = 0, 2 operators, respectively [22][23][24]. The KLOE-2 collaboration, in the most precise measurement of the η → π + π − π 0 Dalitz plot to date, reports all three asymmetry parameters to be consistent with zero [25], superseding many earlier experimental investigations [22,[26][27][28][29][30][31]. Alternatively, C-violation in the phenomenological expansion of the Dalitz-plot distribution, i.e., a two-dimensional Taylor series around its center, can be studied by allowing for both C-conserving and C-violating terms. Until now the KLOE-2 collaboration has probed the first four C-violating terms of this parameterization, which again are all consistent with zero [25]. Thus, experimentally there is no evidence found for C-violation in η → π + π − π 0 .
However, since the 1960s C-violation in this decay has been mostly neglected by theory until recently a new theoretical formalism was proposed in Ref. [21]. In this framework the decay amplitude is decomposed into three contributions that can be associated with operators describing the isospin transitions ∆I = 0, 1, 2. While the Standard-Model contribution is driven almost exclusively by the ∆I = 1 contribution (ignoring isospin breaking of higher order that is known to have only tiny effects [35,36]), the additional BSM amplitudes arise from ∆I = 0, 2 transitions. The individual strengths of the latter are given by two complex-valued normalizations. Physically this approach is more meaningful compared to simple phenomenological (i.e., polynomial) parameterizations, as it allows for a direct extraction of the coupling strengths that may subsequently be matched to underlying BSM operators. The energy dependence of the C-violating amplitudes in Ref. [21] is based on the well-known one-loop representation of the SM decay in chiral perturbation theory (χPT) [37]. The authors find the BSM normalization of the ∆I = 0 amplitude to be between two and four orders of magnitude less rigorously constrained than the ∆I = 2 one, which is a result of the predicted kinematic suppression of the ∆I = 0 transition [34], but again there is no hint for C-violation in η → π + π − π 0 as both BSM normalizations are consistent with zero.
A more rigorous construction of the BSM amplitudes consistent with the fundamental principles of analyticity (a mathematical description of causality) and unitarity (a consequence of probability conservation) can be achieved with techniques from dispersion theory, using the so-called Khuri-Treiman representations [38]. As Sutherland's theorem [39,40], a statement of current algebra, and χPT calculations [35,41] proved that electromagnetic effects are tiny compared to isospin breaking due to the light quark mass difference m u −m d , modern dispersion-theoretical studies of the SM contribution η → 3π [42][43][44][45][46][47][48] focus on a consistent, non-perturbative description of the final-state interactions with the goal to provide information on these fundamental SM parameters. Such a treatment of final-state interactions can also be incorporated in the C-violating amplitudes by establishing the corresponding dispersion relations for the ∆I = 0, 2 transitions. As a by-product, such dispersive amplitude representations allow us to argue more rigorously why the dependence on yet unknown short-distance operators can be subsumed in a single unknown multiplicative constant for each isospin.
The opportunity to investigate C-odd effects as an interference in a Dalitz plot exists similarly for the decay η ′ → ηπ + π − (although without the potential benefit of the SM decay being suppressed by isospin). This is particularly interesting as the possible asymmetry in the distribution of the charged pions in this decay is sensitive to a different class of Cviolating operators from those constrained in η → π + π − π 0 , namely the ones with ∆I = 1. Both decays therefore provide orthogonal probes as far as the isospin structure of the Cviolating operators is concerned. For η ′ → ηπ + π − such an operator must include two derivatives and explicitly reads (1. 2) The experimental limits on the left-right asymmetry A LR and the C-odd contributions of the phenomenological Dalitz-plot expansion measured by the BESIII collaboration [49] vanish again within one standard deviation. Prior to that, measurements by VES [50] as well as an earlier BESIII result [51] came to the same conclusion, albeit with much lower accuracy. While the theoretical description of the SM contribution relying on a sophisticated dispersion-theoretical approach was first established in Refs. [52,53], the incorporation of C-violating effects is missing so far.
In this work we generalize the dispersion-theoretical analysis of C-conserving SM η → 3π and η ′ → ηππ decays to additional C-violating BSM contributions. Accordingly, the presented dispersive representations account for a consistent resummation of the respective three-particle final-state interactions in all allowed isospin transitions. To establish dispersion relations for the new C-violating contributions we split our analysis into two parts, Sect. 2 dealing with η → 3π and Sect. 3 with η ′ → ηππ, and follow the same general strategy in both of them. We start with the definition of the T -matrix elements and the general kinematics of the respective process in Sects. 2.1 and 3.1. In Sects. 2.2 and 3.2 we decompose the amplitudes into ones depending on one Mandelstam variable only, tremendously simplifying the evaluation. These single-variable amplitudes (SVA) are constrained by elastic unitarity as described in Sects. 2.3 and 3.3. Sections 2.5 and 3.5 describe how to extract coupling constants for the effective BSM operators in terms of the subtraction constants. The latter are the free parameters of our dispersive representation, which are fixed by a χ 2 -fit to data in Sects. 2.6 and 3.6. Afterwards we compare our representations of the three-body amplitudes to measurements of the corresponding Dalitz-plot distributions and theoretical constraints in Sects. 2.7 and 3.7. Section 2.8 contains a brief comment on how to generalize the analysis for η → π + π − π 0 to η ′ → π + π − π 0 . We conclude our study with a summary covering both parts in Sect. 4.

Dispersive representation of η → 3π
As our first step towards the analysis of TOPE forces, we will investigate the decay η → π + π − π 0 . For this purpose we rely on the sophisticated and well-established Khuri-Treiman framework [38], in which a set of integral equations for the scattering process ηπ → ππ is established. For the corresponding dispersion relations we only take the dominant elastic pion-pion rescattering into account. With an analytic continuation of the decay mass as well as the Mandelstam variables one can project onto the physical realm of the decay, thereby taking final-state interactions to all orders in perturbation theory into account, and at the same time obtain a manifestly unitary amplitude.
Before going into more detail, let us have a look at the general properties of the η → π + π − π 0 amplitude. Regarding the involved quantum numbers, Bose symmetry demands that C = (−1) I+1 [23], where I is the total isospin of the three-body final state, which has to be distinguished from the isospin of the decaying meson. In analyses consistent with the symmetries of the Standard Model [36,37,[42][43][44][45][46][47][48][54][55][56][57][58][59], the decay amplitude is exclusively driven by isospin-breaking effects. The reason for this lies in the fact that this decay breaks G-parity, whose prerequisite is that either isospin or charge conjugation symmetry is broken, or both. As Standard-Model analyses consider the latter the more cherished symmetry (disregarding the weak interactions), the corresponding amplitudes solely contain ∆I = 1 transitions. 2 On the contrary, in this work we allow for even isospin transitions ∆I = 0, 2 and hence imply C-violation. Moreover, considering that all involved particles are pseudoscalars, the decay at hand preserves parity and one can conclude that CP has to be violated, too. In summary, the C-violating mechanisms are driven by isoscalar ∆I = 0 or isotensor ∆I = 2 operators [23,24,34,60], such that the generalized η → π + π − π 0 amplitude has to be of the form [21] which is split into a contribution for each total isospin denoted by the respective index. In accordance with Refs. [45,47], we factorized out the isospin-breaking normalization of the SM amplitude in terms of the pion decay constant F π and the QCD kaon mass difference. The isoscalar amplitude M C 0 is isospin-conserving but C-violating, the Standard-Model amplitude M C 1 is isospin-violating but C-conserving, and the isotensor contribution M C 2 violates both quantum numbers. Note that isospin symmetry is an accidental (approximate) symmetry of the strong interactions due to the smallness of the two lightest quark masses (as well as their difference) on typical hadronic scales; as we do not know anything about the isospin structure of the BSM operators, there is no reason to assume isospin to be a useful symmetry for them, too, and hence imply any kind of hierarchy between isoscalar and isotensor Cviolation on the underlying, fundamental level.
As a further consequence of Bose symmetry, the C-violating operators can only contribute to the charged decay mode, but not to η → 3π 0 . The latter is thus solely given in terms of M C 1 and explicitly reads as demanded by isospin symmetry. Corrections to Eq. (2.3) arise only due to higher-order corrections such as virtual-photon effects or the charged-to-neutral pion mass difference [35,36].
As the Standard-Model contribution M C 1 has already been extensively studied using Khuri-Treiman equations in Refs. [42-48, 55, 56], in this section we generalize the dispersive analysis by elaborating on the C-and CP -odd amplitudes M C 0 and M C 2 .

Kinematics
Let us define the η → π + π − π 0 transition amplitude in the common manner Up to the overall isospin-breaking normalization, we work in the isospin limit, i.e., M π ≡ M π ± = M π 0 , and conventionally write the corresponding Mandelstam variables as fulfilling the relation Note that the amplitude M C 1 is symmetric under t ↔ u, while M C 0 and M C 2 are both antisymmetric under the exchange of these two Mandelstam variables. In the two-pion center-of-mass system, t and u can be expressed in terms of s and the s-channel scattering angle z s by where σ(s) = 1 − 4M 2 π /s and λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + xz + yz) denotes the Källén function.

Reconstruction theorem
To avoid the intricate analysis of complex functions depending on multiple variables one can exploit a decomposition of the amplitude into single-variable functions. Such a decomposition is commonly referred to as reconstruction theorem. It was first proven that the latter holds exactly up to and including two-loop order in the framework of χPT for ππ scattering [61], followed by generalizations for unequal masses [62] and scattering of mesons belonging to the pseudoscalar octet [63].
Neglecting the discontinuities of D-and higher partial waves, one may express each amplitude of total isospin, on the right-hand side of Eq. (2.1), in terms of functions depending on only one kinematical variable, the relative angular momentum and isospin of the ππ intermediate state [21,55,56]: (2.9) Due to Bose symmetry, the isospin I of the two-pion state fixes the partial wave ℓ unambiguously by means of A 0,2 ≡ A ℓ=0 I=0,2 and A 1 ≡ A ℓ=1 I=1 , with A ∈ {F, G, H}. Note that the single-variable functions F, G, and H are completely decoupled and can be evaluated independently. Furthermore, each single-variable function A I (s) has only a right-hand cut. At this point the charge asymmetry in the C-odd contributions, stemming from the exchange t ↔ u, becomes evident. It is worth noting that the decomposition into singlevariable functions presented here is not unique. The relation between the Mandelstam variables given in Eq. (2.6) allows us to shift the amplitude by polynomials in s, t, and u, i.e., A I → A I + ∆A I , without affecting the reconstruction theorems. The five-parameter ambiguity for the Standard Model reads [47] ∆F 0 (s) = −4a while for the C-odd contributions we find The invariance groups of the amplitudes M C 1 , M C 0 , and M C 2 , given by polynomial ambiguities, are different and independent of each other (in contrast to the erroneous assumption made in Ref. [21]).
Finally, we would like to address the issue of corrections to the reconstruction theorems for M C 1 , M / C 0 , and M / C 2 stated in Eq. (2.9). The next discontinuities, beyond those in Sand P -waves, would come from D-(for even isospin) and F -waves (for odd isospin). Since the symmetry structure of the isoscalar amplitude does not allow for even partial waves, it is obvious that its reconstruction theorem actually holds up to corrections due to F -and higher (odd) partial waves. Moreover, possible D-wave contributions to the discontinuity of the isotensor amplitude are only allowed to have I = 2, which is a nonresonant and extremely small partial wave at low energies. As the validity of the reconstruction theorem for the C-conserving amplitude M C 1 , which neglects discontinuities due to I = 0 D-wave pionpion rescattering, is well-established and tested against very accurate data, we conclude that corrections to the decomposition of the C-violating amplitudes M / C 0 and M / C 2 are necessarily even smaller and therefore entirely negligible.

Elastic unitarity
Within the scope of this work we will exclusively study the dominant elastic rescattering effects, i.e., we restrict the evaluation of the single-variable functions to ππ intermediate states only. In order to obtain an amplitude with manifest unitarity, each single-variable function has to obey the discontinuity relation 3 disc A I (s) = 2i θ(s − 4M 2 π ) A I (s) +Â I (s) sin δ I (s) e −iδ I (s) . (2.12) Here we introduced the so called inhomogeneitiesÂ I (s) that do not have a discontinuity along the right-hand cut and can be evaluated by a projection onto the respective partial wave 4 a I (s) = A I (s) +Â I (s) . (2.13) Note that the full information about the discontinuity of the partial wave along the righthand cut is contained in the respective A I (s). Let us now, for the sake of simplicity, define the angular average (2.14) This allows to write the inhomogeneities for the Standard-Model amplitude in the shortened formF 15) and the ones for the C-violating contributions aŝ (2.16) Note that the argument of A I in Eq. (2.14) is Mandelstam t, meaning that the inhomogeneity in the s-channel single-variable function is determined by contributions of the crossed channels. In other words,Â(s) contains left-hand-cut contributions to the respective partial wave.
In order to obtain a unique solution for the discontinuity relation in Eq. (2.12) it is appealing to first consider the homogeneous case by settingÂ I (s) = 0. 5 The homogeneous 3 Note that for elastic two-body scattering the discontinuity can be replaced by the imaginary part of AI (s). In contrast, for a three-body decay, which is obtained by analytic continuation, the right-hand side of the unitarity equation above is not purely imaginary asÂI (s) becomes complex and Schwarz' reflection principle is not applicable anymore. Therefore we will solely refer to the discontinuity of AI (s). 4 We define the partial-wave projection for a scalar 2 → 2 scattering amplitude TI (s, zs) for fixed twoparticle isospin I by where P ℓ denote the Legendre polynomials. For details on how these partial-wave amplitudes relate to M C 0 , M C 1 , and M C 2 , we refer to Ref. [64]. 5 This scenario is consistent with Watson's final-state theorem, which states that the phase of AI coincides with the phase shift of elastic ππ rescattering. solution corresponds to the one of a pion form factor (of the appropriate quantum numbers) and is given in terms of the Omnès function [65] . (2.17) Using the latter, the general solution becomes where P n−1 (s) is a polynomial in s of order n − 1. Its coefficients are known as subtraction constants, which are the only free parameters of our amplitude. Throughout this paper, we will assume all subtraction constants within the same decay amplitude representation to be relatively real. This is not rigorously true to arbitrary precision, as the SVAs do not fulfill the Schwarz reflection principle, and their discontinuities are complex. However, the potential imaginary parts of the subtraction constants scale with the available three-body phase space, and therefore are tiny for decays such as η → 3π or η ′ → ηππ. This has been tested explicitly for η → 3π [47], making use of the two-loop representation in chiral perturbation theory [59], with the result that imaginary parts in the dispersive subtractions are entirely negligible. This is, however, not true any more for Khuri-Treiman representations of three-body decays with larger energy releases, see φ → 3π [66] or certain D-meson decays [67,68].
Besides these degrees of freedoms, which have to be fixed by data regression or matching to effective theories, the only input for the dispersive η → 3π amplitude are the ππ scattering phase shifts δ I (s).

Subtraction scheme
Choosing the number of subtraction constants n is a rather sensitive issue. Having a purely mathematical look at the dispersion integral in Eq. (2.18), the minimal number is the one at which convergence is ensured. Any additional subtraction just leads to a rearrangement of the equation. The thereby introduced subtraction constants have to fulfill the corresponding sum rule, such that the minimally and higher subtracted integrals are analytically the same. Any deviation from the respective sum rule violates the initially assumed highenergy behavior and is inconsistent as a matter of principle. But allowing the additional subtraction constants to vary from the sum rule suppresses the hardly-constrained highenergy behavior of the dispersion integral and introduces additional degrees of freedom in a fit to experimental data.
Past studies in the same Khuri-Treiman framework as presented here [45,47] put their main focus on maximal precision of the low-energy representation of the η → 3π Standard-Model decay amplitude, and therefore incorporated a rather generous number of subtraction constants. Our aim here is slightly different: we will demonstrate that with rigorous assumptions on the high-energy behavior, and accordingly a minimal number of free parameters, we are still able to describe the Dalitz plot data sufficiently well. Subsequently, we impose the same high-energy asymptotics on the two C-violating amplitudes, and show 0.08 0.1 0.12 0.14 0.16 Figure 1: The S-and P -wave ππ scattering phase shifts δ 0 (red), δ 1 (green), and δ 2 (blue) covering low-energy uncertainty bands as determined by Roy equation analysis [69,70].
Left panel: behavior of the phase shifts in the low-energy region below the KK-threshold at about 50M 2 π . The phase space for η → 3π is indicated by the gray region. Right panel: magnification of the physical decay region.
that as a result, they can be written in terms of one single subtraction constant each. In this manner, we can prove that the mere assumption to describe the BSM amplitudes in terms of a multiplicative normalization only [21] can be justified more rigorously in terms of their analytic properties.
In order to investigate the convergence of the dispersion integral, some assumptions have to be made for the asymptotics of δ I (s) and A I (s). We rely on a Roy equation analysis [69,70] to fix our phase shifts very precisely in the low-energy range, i.e., below about 1GeV 2 , as shown in Fig. 1. Unfortunately the high-energy behavior is not severely restricted by these equations. Therefore we suppose that in the limit s → ∞ the phase shifts approximate constants and analytically continue them accordingly. These limits directly fix the asymptotics of the Omnès functions, which behave for s → ∞ like s −k if δ I (s) → kπ. Further, in order to use the minimal number of subtraction constants, we assume that our amplitudes scale in the limit of large momenta as 20) and are thus even more restrictive than suggested by the Froissart-Martin bound [71]. Finally, with our minimal subtraction scheme we can obtain for the C-conserving Standard- This representation hence depends on three free parameters (all of which are chosen to be real), where Refs. [45,47] employed six. Similarly rigorous schemes with few parameters have previously been suggested in Refs. [43,44,46]. Analogously the C-violating contributions become Note that these representations are not unique. We exploited the ambiguity of the dispersive representation, as given in Eqs. (2.10) and (2.11), to express the single-variable functions in terms of independent subtraction constants only. Conventionally, we shifted the polynomials ∆A I such that the I = 2 amplitudes do not contain subtraction constants. Further, we would like to remark that the normalization of each amplitude of total isospin in Eq. (2.1) has a phase that is fixed unambiguously by T violation and hermiticity. Hence, the subtraction constants ε and ϑ, which absorb these normalizations, are complex quantities with a fixed phase, resulting in a total of five degrees of freedom for M. We furthermore note that, were it not for strong rescattering phases, interference between the (then purely real) Standard-Model and the (purely imaginary) C-violating amplitudes would vanish altogether, and no Dalitz-plot asymmetries would be generated at all. This further increases the importance to implement these rescattering effects via the corresponding phase shifts exactly, using dispersion theory. 6 The fact that these two subtraction constants indeed merely serve as overall normalizations of the BSM contributions becomes evident when realizing that the dispersive representation is linear in the subtraction constants. This very powerful property allows us to write where ν and µ denote generic subtraction constants and M ∈ {M C 0 , M C 1 , M C 2 }. This procedure simplifies the numerical computation tremendously, as the corresponding basis amplitudes A ν I , which obey analogous relations as the M ν in the equation above, can be evaluated once and for all before fixing the subtraction constants.

Taylor invariants
Any interpretation of subtraction constants, which do not have any physical meaning on their own, should be made with caution, as they depend on the chosen subtraction scheme, on the ambiguities of the dispersive representation, and on the not-well-restricted highenergy behavior of the dispersion integrals. Changes in any of the listed aspects are absorbed in the subtraction constants when fitting to data. 7 To overcome these issues we follow the idea of Refs. [45,47], where certain linear combinations of the subtraction constants for the SM contribution were introduced, which are identified as so-called Taylor invariants. To access those, the single-variable amplitudes Inserting the series into the reconstruction theorem for the SM amplitude, cf. Eq. (2.9), one obtains 8 with the Taylor invariants (2.26) These can be used as theory constraints to the SM amplitude when considering that oneloop χPT [37,47] predicts them to be 27) where F 0 was used as an overall normalization by means of f i ≡ F i /F 0 and will furthermore serve to normalize M C 1 . 9 7 Even a simple estimation of the relative and overall size of the two C-odd amplitudes from the subtraction constants ε and ϑ, as in similar fashion assumed by Ref. [21], may be misleading. Notice that an apparent difference in these coefficients can be due to the compensation of the relative, arbitrary normalization of the basis solutions G ε I and H ϑ I when comparing to data. 8 For simplicity, here and in the following we denote the order of neglected higher-order polynomial terms by O(p 2n ), which should not be confused with the counting scheme of the chiral expansion that may include nonanalytic dependencies on quark masses etc. 9 In principle one can also define Taylor invariants for the SM amplitude at the two-loop level in χPT [48,59]. However, as demonstrated in the analysis of Ref. [47], a high-precision matching requires a more flexible dispersive amplitude (i.e., more than three subtraction constants). Aside of this small flaw, we will demonstrate that our dispersive representation of the SM amplitude describes both the experimental Dalitz-plot distribution and the one-loop chiral constraints very well.
We now apply the same strategy to the C-odd contributions. The effective BSM operators of Eq. (1.1), which arise from elementary considerations such as crossing symmetry and the correct behavior under time reversal, demand the amplitudes for the ∆I = 0 and ∆I = 2 transitions at lowest contributing order to be of the form where the couplings have the dimensions [g 0 ] = GeV −6 and [g 2 ] = GeV −2 , respectively. It has to be remarked that this simple polynomial expansion is by far less accurate than the full dispersive representation, but allows one to match the couplings in a convenient way.
Reproducing the structure in Eq. (2.28) with the Taylor series from above we obtain which we wrote in explicit dependence on the subtraction constants using the Taylor invariants for the basis amplitudes A ν I , by means of C G 1 = ε C G ε 1 etc. In this form, T violation demands the coupling constants g 0 and g 2 to be real-valued. To satisfy this condition, the subtraction constants must be proportional to the complex conjugate of the linear combinations of Taylor invariants, by means of as an example for ε. While this condition fixes the phase of ε, the constant c ε is left as the only degree of freedom. We proceed analogously with ϑ. As we extract the Taylor invariants by an expansion of the single-variable amplitudes around s = 0, we are well below the dipion threshold, such that the contributions of the dispersion integrals are negligible for the decay at hand. Consequently, we can drop the real part of the subtraction constants, which have no visible effects on observables.

Fixing the subtraction constants
Once the basis solutions A ν I for the Khuri-Treiman coupled integral equations are evaluated numerically, one can determine the free parameters of our dispersive representation for η → 3π. In summary we have the subtraction constants α, β, γ for the SM amplitude, where one of these can be seen as an overall normalization, as well as ε fixing the C-violating isoscalar contribution and ϑ for the isotensor one.
To determine these degrees of freedom we employ a χ 2 -regression to three different data sets: • the Dalitz-plot distribution of η → π + π − π 0 from the KLOE-2 collaboration [25], • the Dalitz-plot distribution of η → 3π 0 from the A2 collaboration [72], and • the Taylor invariants of M C 1 from one-loop χPT.
Note that the latter two only address the three free parameters of the SM amplitude. However, these two data sets help to fix the relative phases between the contributions of different total isospins, and furthermore any shift in M C 1 may affect the BSM contributions when additionally comparing to the data of Ref. [25].
Let us first turn our attention to the experimental data sets from the KLOE-2 and A2 collaborations. The KLOE-2 collaboration provides the world's highest statistics for the measurement of the Dalitz-plot distribution in η → π + π − π 0 . The data distributes about 4.6 · 10 6 events over 371 bins, where all bins overlapping with the physical boundaries were discarded. On the other hand, the A2 collaboration provides altogether 441 bins for a single Dalitz-plot sextant, exploiting the symmetry of η → 3π 0 , and accepts bins overlapping with the phase space boundary. It supersedes many earlier experiments on η → 3π 0 , which mostly concentrated on the leading nontrivial Dalitz-plot slope parameter [31,[73][74][75][76][77]. Let us refer to the experimental Dalitz-plot distributions by D exp c,n , where the index denotes the charged or neutral channel, respectively. The binning is given, as commonly done, in terms of the dimensionless and symmetrized coordinates x i c,n , y i c,n , where the additional index denotes the i-th bin at its center. These explicitly read The indices labeling the Mandelstam variables correspond to the respective kinematic map given in Ref. [47]. We compare the experimental measurements to the dispersive Dalitz-plot distributions by integrating our amplitudes M c,n over the respective bin The discrepancy functions χ 2 c,n for the charged and neutral data sets are then defined by (2.33) To build in theory constraints on the Taylor invariants of the SM amplitude from one-loop χPT we introduce [47] where the f χPT i denote the theoretical predictions listed in Eq. (2.27). To define a realvalued discrepancy function we restrict our analysis to the real parts of the f i and discuss the effects of their imaginary parts in Sect. 2.7.1.
When carrying out the combined regression to all three data sets we minimize the combined discrepancy function  Table 1: Summary of the four considered fit scenarios: SM c (exclusive, C-conserving), BSM c (exclusive, C-violating), SM tot (combined, C-conserving), BSM tot (combined, Cviolating). For fits obtained by dropping the contributions of χ 2 0 and χ 2 n to the total discrepancy function χ 2 tot , their values are put in brackets. All values refer to fit results of our central solution.
and fix the normalization of M c such that it reproduces the Taylor invariant F 0 from Eq. (2.27). Before fixing the subtraction constants, one has to consider higher-order isospin corrections due to the mass difference of the neutral and charged pions in the final state to obtain an accurate description of the experimental measurements. To this end, we follow the same strategy as proposed in Ref. [47], which shall serve as a reference for explicit formulas. The dominant isospin-breaking contribution can be taken into account by a kinematic map, such that the boundaries of the Dalitz plot in the isospin limit are mapped to the ones for physical masses. All remaining isospin-breaking effects are assumed to be mostly absorbed by electromagnetic correction factors K c,n for the charged and neutral decay modes of η → 3π resulting from one-loop representations in χPT [35]. While the kinematic map will be applied to both the C-even and C-odd amplitudes, K c,n only enter the SM amplitudes, as we are yet missing any effective theory to account for analogous corrections in the Cviolating amplitudes. Due to the absence of I = 0 S-wave contributions, in general we expect such electromagnetic effects to be even smaller in that case.
When minimizing the χ 2 as described above, we distinguish between the four scenarios • SM c : exclusively minimize χ 2 c with M C 0,2 = 0 , • BSM c : exclusively minimize χ 2 c with the full amplitude M c , • SM tot : minimize χ 2 tot with M C 0,2 = 0 , • BSM tot : minimize χ 2 tot with the full amplitude M c .
A summary of the individual χ 2 contributions to the four scenarios is given in Table 1.
We find for all fit scenarios considered a good agreement of our dispersive amplitude with data. Overall the individual parts of the discrepancy function χ 2 0 , χ 2 c , and χ 2 n in the four different scenarios are almost identical. In fact, the dispersive representation is already perfectly fixed by the KLOE-2 data on η → π + π − π 0 alone, with the η → 3π 0 Dalitz-plot distribution and the Taylor invariants for M C 1 being a prediction in excellent agreement with  data. Accordingly, the differences between the results for the exclusive and combined fits are marginal, i.e., comparing SM c vs. SM tot and BSM c vs. BSM tot . Taking the C-violating contributions into account and comparing SM c vs. BSM c or SM tot vs. BSM tot respectively, we find a minor improvement of χ 2 c by about 1.1%, whereas χ 2 0 and χ 2 n do not change at all at the given level of accuracy. Furthermore, comparing the resulting discrepancy functions of the KLOE-2 and A2 data sets, we notice a slightly worse description of the Dalitz-plot distribution for the neutral η → 3π 0 mode. This small tension of the dispersive representation for M n and the experimental measurement from A2 has also been observed in Ref. [47]. Nevertheless, the experimental data of both the charged and neutral mode together are well described. Consequently, adding the contributions of M / C 0 and M / C 2 to our dispersive representation for M c has no visible effect on the determination of M C 1 . Let us now elaborate on the error analysis. For the latter we consider the experimental uncertainties from the KLOE-2 and A2 Dalitz-plot distributions, 10 the uncertainty originating from χPT constraints including the Taylor invariants for M C 1 (2.27) and the electromagnetic correction factors K c,n from Ref. [47], and the uncertainty resulting from the variation of the phase shift input in the low-and high-energy region, cf. Fig. 1. We will treat all these sources of error as symmetric and Gaussian distributed. Accordingly, the combined total uncertainties are found by adding the individual contributions in quadrature and the presented correlation matrices are calculated from the respective total covariance matrices of the investigated quantities.
For the sake of completeness we quote the subtraction constants determined with the combined regressions SM tot and BSM tot in Table 2, which underline the findings pointed out in the previous paragraphs, and illustrate the corresponding comparison to the KLOE-2 Dalitz plot in Fig. 2. Due to the reasons stated in Sect. 2.5 we will refrain from any further discussion of the subtraction constants and instead have a look at actual observables in the following sections. Based on the observations discussed above, we will henceforth exclusively refer to the results obtained with the scenario BSM tot .
After fixing the C-conserving contribution with the subtraction constants listed in Table 2 we can compare the three contributions M   3. With the shown extrapolation to the region of the ρ(770) resonance, we observe that the isovector and isotensor contributions, i.e., F 1 and H 1 , are in good agreement with each other over an energy range exceeding the physical decay region. On the contrary, these SVAs are significantly different from the isoscalar one, i.e., G 1 . 11 Hence the approximation that F I (s) = G I (s) = H I (s) as assumed in Ref. [21] may not be accurate enough to investigate possible future measurements of C-violating effects in η → π + π − π 0 .

Extraction of observables
In this section we present the numerical results for several observables that can be extracted from our dispersive representation of the η → 3π amplitudes. We start our discussion by first investigating theoretical and experimental constraints imposed on the SM amplitude and focus on a comparison of our results to the established analysis of Ref. [47], which shall serve as a consistency check of our dispersive representation. We show that our minimal subtraction scheme for the SM amplitude meets these requirements and can thus argue that the application of this subtraction scheme to the BSM amplitude, cf. Eq. (2.22), is justified.
Subsequently we have a closer look at C-violating observables of the η → π + π − π 0 Dalitz-plot distribution, the occurring asymmetries, and the coupling strength of effective BSM operators with isospins ∆I = 0 and ∆I = 2.  Table 2. The phase space for η → 3π is indicated by the gray region.

Standard Model constraints
Let us start the discussion concerning the validity of our SM amplitude M C 1 by having a look at theoretical constraints from one-loop χPT. For this purpose we extract the Taylor invariants as described in Sect. 2.5. Our value for the normalization of the Taylor invariants yields F 0 = 1.176(53) − 0.0094 (14) i. As stated previously we fixed the normalization of M C 1 so that it reproduces Re F 0 from Eq. (2.27), but allowed Im F 0 to vary. Since the latter is exclusively generated by contributions of the dispersion integrals to Eq. (2.21) it is roughly two orders of magnitude smaller than Re F 0 and will be neglected from now on. For the real parts of the reduced coefficients f i and their correlation we hence find In contrast to the dispersive representation of Ref. [47], which uses a subtraction scheme for M C 1 involving six independent subtraction constants, our minimalist scheme (2.21) is extremely stiff. Therefore it does not allow for a large variation of the reduced Taylor invariants, cf. Table 1. Similar to Im F 0 the imaginary parts of the reduced invariants Im f 1 = 0.193(29) GeV −2 , Im f 2 = −0.006(85) GeV −4 , and Im f 3 = −0.128(39) GeV −4 are found to be small. Next, we want to consider the behavior of M C 1 at its soft-pion point, i.e., in the limit where the four-momentum of one of the pions vanishes. As current algebra dictates, the amplitude M C 1 must exhibit a zero at this point. In terms of the Mandelstam variables we will find two of these so called Adler zeros at s A = t A = 0 and s A = u A = 0 which is in perfect agreement with the χPT prediction. Nevertheless, we want to mention that the Adler zeros are shifted slightly away from the critical lines s = t and s = u. The dominating error source in Eq. (2.37) stems from the low-and high-energy uncertainties of the phase shift input, cf. Fig. 1. One last consistency check regards the observables of the neutral channel η → 3π 0 . Due to the symmetry under exchange of any two Mandelstam variables we stick to the common phenomenological parameterization in terms of the polar coordinates z n and φ n given by |M n (z n , φ n )| 2 ∼ 1 + 2α z n + 2β z 3/2 n sin 3φ n + . . . . where the slope α agrees well with the Particle Data Group (PDG) world average [78] and the parameter β is compatible with the findings of the A2 collaboration [72] as well as with the dispersive analysis of Ref. [47]. The extraction of higher parameters is beyond the scope of this work. Finally, we can also calculate the ratio BR(η → 3π 0 )/BR(η → π + π − π 0 ), which can be computed from partial decay widths Γ c,n defined by where S c = 1 and S n = 6 denoting the symmetry factors and D c,n the integrals of the Dalitz-plot distributions over the full phase space. Since contributions antisymmetric under t ↔ u cancel, D c is determined entirely by |M C 1 | 2 up to corrections quadratic in the BSM couplings. We extract in perfect agreement with the PDG world average [78]. Note that the uncertainty quoted in Eq. (2.41) is totally dominated by the errors on the electromagnetic correction factor K c,n from Ref. [47]. As our minimal subtraction scheme meets all the presented constraints imposed on the SM amplitude, we conclude that there is no objection when applying it to the BSM contributions.

Dalitz-plot distributions
We now turn our focus to the determination of C-violating observables in the η → π + π − π 0 Dalitz-plot distribution. As already observed in Sect. 2.6, patterns arising from TOPE forces have a vanishingly small influence on the goodness of the regression. Nevertheless, we show to which order of magnitude C-and CP -violating signals in η → π + π − π 0 , as predicted by our dispersive representation, can be restricted with the currently most precise measurement of the respective Dalitz plot [25]. For this purpose, it may be advantageous to decompose the Dalitz-plot distribution of the total amplitude in Eq. (2.1) into its constituents by means of where we neglected all contributions that are quadratic in C-violating amplitudes, i.e., 2 are of the same order of magnitude. Like the SM contribution |M C 1 | 2 , all effects quadratic in C-violation are symmetric under t ↔ u and will therefore not contribute to the mirror symmetry breaking.
Due to the small phase space of the decay at hand the Dalitz plot is typically parameterized by a polynomial expansion around its center, by means of where the coefficients a, b, etc., are called Dalitz-plot parameters. By now, the first seven coefficients of this phenomenological parameterization have been studied by the KLOE-2 collaboration [25]. Note that non-vanishing values of the coefficients c, e, h, and l odd in x c would directly implicate C-violation in η → π + π − π 0 decays. We first access the Dalitz-plot parameters for the C-conserving contribution, generated exclusively by M C 1 , by employing a two-dimensional Taylor  The uncertainties of these four parameters are dominated by the statistical error of the KLOE-2 data, while all other sources of uncertainty do not yield any significant contribution to the error budget. 12 Accordingly, we can confirm that all C-violating parameters vanish within 2σ at most. Furthermore, the C-violating parameters turn out to be at least one order of magnitude smaller than d and g, which are the smallest coefficients of the C-  A comparison of the extracted Dalitz-plot parameters with the results from KLOE-2 as well as the two most recent dispersive analyses on C-conserving η → 3π decays [46,47] are summarized in Table 3.

Asymmetries and BSM couplings
Besides these coefficients, we can also investigate three asymmetry parameters to quantify C-violating effects in the η → π + π − π 0 Dalitz-plot distribution: the left-right A LR , the quadrant A Q , and sextant A S asymmetry parameters [22][23][24]. These asymmetries compare the population of the Dalitz-plot distribution in the different sectors defined by the Dalitzplot geometry, cf. Fig. 5. To quantify these asymmetries we follow Ref. [21] by defining 48) with N = N R + N L and denoting the normalized number of events for the total amplitude within each region C. where all three asymmetry parameters are given in units of 10 −4 . We find A LR , A Q , and A S in good agreement with the results reported by the KLOE-2 collaboration [25]. Again, there is no hint for C-violation as all three asymmetries are compatible with zero in not more than 1.8σ. Note that the error budget in Eq. (2.50) is completely dominated by the statistical uncertainties of the KLOE-2 data. 13 In contrast to experimental studies of C-violating effects in the η → π + π − π 0 Dalitz-plot distribution, which are limited to the investigation of x c -odd coefficients of the phenomenological parameterization (2.43)  Note that for the central values we find a ratio of |g 0 /g 2 | ≈ 10 3 GeV −4 . This can be understood as follows. Generically, as we have remarked above, the operator X / C 0 is kinematically suppressed compared to X / C 2 by 4 orders in the chiral expansion; this means that we would expect their coefficients to behave as |g 0 /g 2 | ∼ (1 GeV) −4 , the scale given by the chiral symmetry breaking scale 4πF π ≈ 1.16 GeV. As the momenta throughout the η → 3π Dalitz plot are of order M π (note the available phase space M η − 3M π ≈ M π ), this would lead to a relative suppression of the isoscalar transition with respect to the isotensor one of roughly (M π /1 GeV) 4 ≈ 4 × 10 −4 . In fact, however, the data constrains both amplitudes including their respective coupling constants about equally, cf. Fig. 5, which means that the experimental sensitivities rather behave like |g 0 /g 2 | ∼ M −4 π ≈ 2.6 × 10 3 GeV −4 , in good agreement with what we observe. This behavior of the amplitudes M / C 0 and M / C 2 has also been observed in Ref. [21].
Furthermore we can utilize these coupling strengths to obtain a more general representation of the Dalitz-plot asymmetries. Carrying out the phase space integrals individually for contributions involving interference effects of M / C 0 or M / C 2 in the Dalitz-plot distribution, we find that the asymmetry parameters (2.50) given in units of 10 −4 are related to the BSM couplings g 0 and g 2 by
In these relations g 0 and g 2 enter in units of 1 GeV Once more, all asymmetry parameters are given in units of 10 −4 .
To conclude the discussion of η → 3π, we would like to comment on the future experimental focus to set more severe bounds on C-and CP -violation. We disrecommend using the polynomial parameterization of the Dalitz plot from Eq. (2.43), which is too inaccurate for this purpose, mostly because the order of the polynomial, i.e., the number of degrees of freedom, is not known a priori and depends strongly on the precision of the measurement. On the other hand, the measurement of two out of the three Dalitz-plot asymmetries is in principle sufficient to fix the two degrees of freedom in our amplitude representation of C and CP violation. Note however that we predict strong correlations between the three asymmetries, which would become even more significant if, as naturalness suggests, the isoscalar contribution is strongly suppressed compared to the isotensor one therein. We therefore advocate the use of the more physical decay amplitudes with proper phase behavior in future experimental analyses.

Generalization to η ′ → 3π
As η and η ′ have largely the same quantum numbers and differ mainly due to their masses, the fundamental decay mechanisms into the 3π final states are also identical. In the Standard Model, η ′ → 3π is also almost exclusively due to the light-quark-mass difference, and the classification in terms of isospin amplitudes works in exactly the same way as for η → 3π. Consequently, the same goes for C-violating decay mechanisms. A major difference concerns only the total widths of η and η ′ : while the partial widths of both mesons into three pions are of comparable size, the lifetime of the η ′ is shorter by about a factor of 150, and hence the branching ratios make η ′ → 3π relatively rare decay modes. As a result, high-precision investigations of the corresponding Dalitz plots on the same level as for η → 3π, with the goal to put limits on C-odd effects therein, will most likely remain extremely difficult in the near future. To date, the BESIII collaboration has investigated the decay dynamics in η ′ → 3π most precisely, with a determination of the respective branching ratios [79], a measurement of the η ′ → 3π 0 Dalitz plot [31], and the first amplitude analysis for both charged and neutral final states [80].
Here, we merely intend to estimate the relative size between isoscalar and isotensor C-violating transitions in η ′ → π + π − π 0 : due to the significantly larger available phase space, we suspect the strong kinematic or chiral suppression of the isoscalar amplitude in η → π + π − π 0 to be lifted to a certain extent; an expectation that will be borne out below. As a result, despite the experimental difficulty due to the smaller branching ratio, as a matter of principle η ′ → π + π − π 0 will be much more sensitive to the isoscalar C-odd operators. For the purpose of this qualitative investigation, it is sufficient to consider a dispersive representation of the η ′ → π + π − π 0 decay amplitude as a rescaled version of η → π + π − π 0 , with the mass of the η replaced by the one for its heavier version η ′ , hence increasing the available phase space. We omit the incorporation of any inelasticities, like via the dominant decay channel η ′ → ηππ.
For the purposes of our rather qualitative argument, we only investigate the phase space distributions of the C-odd contributions, because both amplitudes M C 0 and M C 2 only depend on one complex subtraction constant each. Since both amplitudes are driven by the same type of C-and CP -violating operators as the corresponding ones in η → π + π − π 0 , we suppose for our qualitative estimation that the respective coupling constants g 0 and g 2 are equal in both decays. Under this assumption we adjust the normalization of M C 0 and M C 2 in η ′ → π + π − π 0 in terms of the subtraction constants ε and ϑ to g 0 and g 2 as extracted from the central results of the BSM couplings in Eq. (2.51). Note that the contribution of the dispersion integral to the real part of the subtraction constants is in this case not negligible due to the increased phase space. Hence we fix the subtraction constants according to Eq. (2.30). The thereby generated distributions of the real and imaginary parts of the C-odd amplitudes in Fig. 6 show that the chiral suppression of the isoscalar transition with respect to the isotensor one is attenuated significantly by the increased phase space, such that M C 0 dominates M C 2 by roughly two orders of magnitude. More precisely, we predict that the relative sensitivity to the isoscalar transition is increased by about two orders of magnitude in comparison to the analogous η decay. This scaling can be qualitatively understood: we have emphasized that the isoscalar C-odd operators are suppressed by four orders in the chiral expansion with respect to the isotensor ones. However, the η ′ decay into three pions is far less a low-energy decay: the available phase space is larger by about a factor of (M η ′ − 3M π )/(M η − 3M π ) ≈ 4. Taking this to the fourth power correctly predicts an increased relative sensitivity by roughly a factor of 250.
A more rigorous analysis of C-odd effects can be performed once a dispersion-theoretical fit to η ′ → 3π Dalitz plots within the Standard Model is accomplished [81].

Dispersive representation of η ′ → ηππ
In this section we turn our attention to another class of TOPE forces by studying the decay η ′ → ηπ + π − . Considering the quantum numbers of the involved mesons, one can argue in a similar manner as previously in Sect. 2: as the decay at hand preserves G-parity, transitions of even isospin ∆I = 0, 2 conserve C, while odd ones violate the latter. Thus we can write the most general amplitude up to linear order in isospin breaking as where for this decay, as opposed to η (′) → 3π, the isoscalar amplitude M C 0 is isospinand C-conserving, whereas the M C 1 violates both quantum numbers. Note that the decay Im M / C 2 × 10 3 Figure 6: Estimation for the distribution of the real (left) and imaginary (right) parts of M C 0 (top) and M C 2 (bottom) in η ′ → π + π − π 0 over the allowed phase space. As a rough estimation, the normalizations of each amplitude are fixed by the central fit result obtained for η → π + π − π 0 in Table 2. In contrast to Fig. 5 additional zero lines occur as a consequence of the interference with the ρ-resonance, which lies in the kinematically accessible region. η ′ → ηπ + π − is sensitive to a different class of C-and CP -violating operators from those tested in η (′) → π + π − π 0 , namely the ones for transitions with ∆I = 1.
For the evaluation of the overall amplitudes we again rely on the Khuri-Treiman framework, which was already applied to the Standard-Model contribution M C 0 in Ref. [53]. The set of dispersion relations is built from the two scattering processes η ′ η → ππ (s-channel) and η ′ π → ηπ (t-channel). Once more, we allow only for elastic rescattering. In order to determine the C-odd amplitude we follow the same agenda as laid out in Sect. 2.

Kinematics
Define the η ′ → ηπ + π − transition amplitude as usual by For the invariant masses we stick to the convention These Mandelstam variables satisfy the relation For the s-channel scattering amplitude η ′ η → ππ one may write For the t-channel η ′ π → ηπ we have Using the kinematic function we can express the t-channel scattering angle as As a consequence of crossing symmetry the corresponding relations for the u-channel can be obtained by exchanging the variables t ↔ u and z t ↔ −z u . The scattering channels have the physical thresholds (3.10)

Reconstruction theorem
In the ongoing, we restrict our amplitude to discontinuities in the lowest contributing partial waves, i.e., to ℓ = 0 for ππ states with isospin I = 0, or ℓ = 1 for those with I = 1, and to ℓ = 0 for the ηπ system with I = 1. We neglect the phase of the ηπ P -wave, which has exotic quantum numbers (i.e., no resonances are expected in the quark model), and is as suppressed at low energies in the chiral expansion as D-and higher partial waves [82]. With these approximations the decomposition of the Standard-Model amplitude in terms of single-variable functions takes the simple form [53] M C 0 (s, t, u) = F ππ (s) + F ηπ (t) + F ηπ (u) , (3.11) with the abbreviations F ππ (s) ≡ F ℓ=0 I=0 ππ (s) and F ηπ (t) ≡ F ℓ=0 I=1 ηπ (t). In this notation the indices ππ and ηπ denote the two-particle final state of the respective scattering process. In a similar fashion we obtain the reconstruction theorem for the C-violating amplitude which can be read off along the lines of Ref. [64]. In this equation we use the short form G ππ (s) ≡ G 1 1 ππ (s) and G ηπ (s) ≡ G 0 1 ηπ (s). The ambiguities of these representations are given by the transformations which leave the full amplitudes unaffected.

Elastic unitarity
To ensure the conservation of probability, the single-variable functions have to obey (3.14) with A ∈ {F, G} and the indices of the phase shifts labeling the respective two-particle intermediate states. Note that in case of M 0 the ππ-state has isospin I = 0, such that δ ππ = δ I=0 ππ for A = F. Analogously, the C-odd contribution M C 1 is driven by a ππ-state that has isospin I = 1, i.e., δ ππ = δ I=1 ππ for A = G. Introducing the abbreviations the inhomogeneities for the Standard-Model amplitude, obtained by a partial-wave projection as described in Sect. 2.3, becomê 16) and the ones entering the C-violating amplitude yield (3.17) Analogously to Eq. (2.18) we can write the general solutions as with two distinct subtraction polynomials P n−1 ππ and P n−1 ηπ of order n − 1. The index of each Omnès function decides which scattering phase shift is used according to Eq. (2.17). In addition to that, one has to differentiate the case Ω ππ = Ω I=0 ππ for A = F from Ω ππ = Ω I=1 ππ for A = G. As our numerical input, we use the same ππ phase shifts as detailed in the discussion of η → 3π in Sect. 2.4. For the ηπ S-wave, we employ the phase of the corresponding scalar form factor constructed in Ref. [83], further refined by imposing constraints from γγ → ηπ 0 [84]. This phase, including the associated uncertainties, is shown in Fig. 7.

Subtraction scheme
In this section we proceed in the same fashion as in Sect. 2.4 to fix the yet undetermined number of subtractions entering the dispersive representation in Eq. (3.18). We assume that the involved phase shifts behave in the limits s → ∞ or t → ∞, respectively, as Furthermore, we demand the asymptotics of the single-variable functions F ππ and F ηπ resulting from the Froissart-Martin bound [71] F ππ (s) = O(s) , This results in a representation of the corresponding SVAs involving four (real) subtraction constants, (3.21) In Ref. [53], a more rigorous scheme with asymptotics analogous to those discussed for η → 3π in Sect. 2.4 and, correspondingly, less subtractions was employed in parallel, and found to describe the Dalitz plot data similarly well, while being more susceptible to sizeable uncertainties due to high-energy input to the dispersion integrals. With the adjusted input for ηπ scattering [84], this reduced scheme ceases to work well [81]. We regard this partly as an artifact of the extremely slow asymptotic rise of the ηπ phase shift, cf. Fig. 7, and therefore decide to stick to the more restrictive asymptotics for the C-odd contribution all the same, in order to avoid a proliferation of subtraction constants therein. The assumptions for G ππ and G ηπ hence are such that the resulting C-violating SVAs are given by (3.23) Conventionally, the polynomial ambiguities from Eq. (3.13) were shifted such that a minimal number of subtraction constants contributes to the A ηπ . Again, the phase of the subtraction constants ̺ and ζ is fixed by T violation, so that M C 1 has two real-valued degrees of freedom, in contrast to the C-violating isoscalar and isotensor contributions in η → 3π which are fixed by a single normalization each. The numerical implementation proceeds in analogy to the strategy presented in Sect. 2.4.

Taylor invariants
As pointed out in Sect. 2.5, the subtraction constants fixing our dispersive representation are no suitable observables. Therefore we again introduce their linear combinations as ambiguity-free Taylor invariants obtained by an expansion of the SVAs around s, t = 0, i.e.,

(3.24)
Of course the series coefficients take different values for SM and BSM contributions. Applying these expansions to the reconstruction theorem (3.11) allows us to express the SM amplitude by where we dropped terms of cubic order in the Mandelstam variables and higher. The BSM operator driving the ∆I = 1 transition as introduced in Eq. (1.2) demands that the matrix element takes the form 27) where in addition to the effective isovector coupling g 1 , we also consider the leading sdependent correction δg 1 . In terms of the Taylor coefficients these quantities read Note that the additional parameter δg 1 ensures that the degrees of freedom of the Taylor expansion match the ones of the dispersive representation for M C 1 . Both couplings are realvalued as demanded by T violation and give rise to the phases of the subtraction constants ̺ and ζ. The latter can be considered as purely imaginary due to the small available phase space.

Fixing the subtraction constants
According to the subtraction scheme chosen in Sect. 3.4, the dispersive representation of the SM amplitude M C 0 contains the four degrees of freedom α, β, γ, λ, where again one subtraction constant can be chosen to fix the overall normalization. The C-violating isovector contribution M C 1 has a total of three parameters ̺, ζ, and ϕ, where the latter fixes the complex phase between M C 0 and M C 1 . After solving for the basis solutions of the dispersive representation in Eqs. (3.21) and (3.23), these subtraction constants can be determined by a comparison to data.
In contrast to Sect. 2.6, we only consider one single data set, i.e., the Dalitz-plot distribution D of η ′ → ηπ + π − from the BESIII collaboration [49]. The latter provides the currently most precise measurement including 3.51 × 10 5 events extracted from J/ψ decays in terms of the symmetrized coordinates with Q η ′ = M η ′ −M η −2M π . We refrain from including data sets on η ′ → ηπ 0 π 0 [49,85,86] in the analysis, as, in contrast to the case of η → 3π 0 , they do not provide truly independent information on the SM amplitude, but rather probe subtle isospin-breaking effects [53,87]. We determine the subtraction constants by minimizing the discrepancy function for which we compute our dispersive amplitude M on the discrete grid covering the centers of all measured bins and normalize M to reproduce the according experimental decay width Γ(η ′ → ηπ + π − ) = 79.9(2.7) keV taken from the PDG [78].
We proceed by carrying out the regression using the pure SM amplitude M C 0 as well as the one for the full BSM contribution M = M C 0 + M C 1 . The results for these fit scenarios, denoted as FIT SM and FIT BSM , are listed in Table 4 and the corresponding subtraction constants can be found in Table 5. We observe that the additional inclusion of the Cviolating ∆I = 1 transition does not have any visible influence on the overall goodness of the regression. As an illustration of the latter we show the phase space corrected x-and  Figure 8: Dalitz-plot projections in x-and y-direction, which are divided by the corresponding phase space. We show the measurement of Ref. [49] overlayed with our dispersive representations covered by the red error bands. Note that the theoretical x-projection on the left is not perfectly symmetric, due to C-violating contributions. In both panels we depict our central solution for the C-conserving part |M C 0 | 2 by the dotted blue line.  y-projections of the Dalitz plot in Fig. 8. Note that the small effects of mirror symmetry breaking are apparent in the x-projection.
Due to the fact that the current constraints for the η ′ → ηπ + π − SM amplitude are by far less restrictive than the ones pointed out for η → 3π in Sect. 2.7.1, we omit an elaborate analysis of the asymmetric systematical errors when varying the input for the ηπ phase shift shown in Fig. 7 [81]. However, we remark that these systematical errors for the SM amplitude may increase up to the same order of magnitude as the corresponding statistical ones. In either way, the C-violating observables in the central scope of this analysis are dominated by their statistical uncertainties.

Extraction of observables
In this section we work out the numerical results of our dispersive representation for various C-violating observables in the η ′ → ηπ + π − amplitudes. Similar to Sect. 2.7 we first discuss the validity of our SM amplitude. To this end we extract the Adler zeros and the Taylor invariants in Sect. 3.7.1. Thereafter we extract patterns of C-violation in the Dalitz-plot distribution, investigate the occurring asymmetries, and compute the coupling strength of an effective isovector BSM operator X C 1 .

Standard Model constraints
The Taylor invariants F i defined in Sect. 3.5 allow us to extract coefficients that can be compared to theoretical analyses for the η ′ → ηπ + π − SM contribution as for instance large-N c χPT or RχT [88,89].
As described in Sect. 3.4, we use four real-valued subtraction constants to fix the SM amplitude. These can be translated to the Taylor invariants where F 0 serves as an overall normalization by means of f i ≡ F i /F 0 . Possible imaginary parts of the Taylor invariants are exclusively generated by the dispersion integrals (3.21) and are disregarded in the following. Furthermore we want to study the behavior of the SM amplitude outside the physical region at its soft-pion points. Chiral SU(2) R ×SU(2) L symmetry expects two Adler zeros to show up at (t−u) = ±(M 2 η ′ −M 2 η ) along the line s = 0 in the limit of massless pions [90][91][92]. Therefore, in analogy to Ref. [53] we study our dispersive amplitude for on-shell pions along the critical line s = 2M 2 π and find two zeros at An updated analysis of the SM η ′ → ηππ decay presented in Ref. [53] is currently in progress [81], based on the latest high-statistics Dalitz-plot measurements from A2 [86] and BESIII [49] for the charged and neutral decay modes.

Dalitz-plot distribution
Let us continue our discussion on C-violating patterns arising from the ∆I = 1 transition η ′ → ηπ + π − Dalitz-plot distribution and quantify corresponding observables. Dropping the dependencies on the coordinates x and y and neglecting the contribution of |M C 1 | 2 , the Dalitz-plot distribution arising from Eq. (3.1) can be written as which is depicted in Fig. 9. We observe a similar, however slightly flattened, hierarchy as in the case of η → 3π worked out in Sect. 2.7.2. The interference term giving rise to the Dalitz-plot asymmetry is constrained to be three orders of magnitude smaller than the SM contribution |M C 0 | 2 , whereas the pure ∆I = 1 contribution |M C 1 | 2 is suppressed by four orders of magnitude. We conclude that the current state of precision for the η ′ → ηπ + π − Dalitz plot merely restricts the effects of the C-violating isovector transition to the relative per mille level.
Given the small phase space of the process, the momentum distribution is quite smooth and commonly approximated by the same expansion as introduced in Eq. (2.43), but with adapted coordinates x and y from Eq. (3.29). The BESIII collaboration finds that the first three C-even coefficients a, b, and d of this expansion are sufficient to parameterize the Dalitz plot, as all other parameters of higher orders in x and y, as well as all parameters odd in x indicating C-violation, are found to be compatible with zero within less than one standard deviation. A two dimensional Taylor  where we neglect correlations smaller than 1% on the right-hand side. Considering the respective errors we find a perfect agreement of all our parameters with the experiment [49].
In particular there is no indication for C-violation as c and e are effectively zero.

Asymmetry and BSM coupling
To finalize our analysis we quantify the asymmetry and the coupling strength of the ∆I = 1 transition in η ′ → ηπ + π − and apply the same procedure as in Sect. 2.7.3. We find the leftright asymmetry in units of 10 −3 to be A LR = 2.1(1.5).

Summary
In this study, we have put the pioneering work of Ref. [21] for C-and CP -violating amplitude representations in the decay η → π + π − π 0 into a rigorous dispersion theoretical framework, and extended the formalism to the analysis of C-and CP -violation in the hadronic threebody decays of the η ′ . Strictly relying on the fundamental principles of analyticity and unitarity, we constructed all three η → π + π − π 0 amplitudes of distinct total isospin, i.e., the SM amplitude M C 1 as well as the C-violating isoscalar and isotensor contributions M C 0 and M C 2 , non-perturbatively based on ππ phase shifts. We demonstrated that the same constraints-all amplitudes are not allowed to grow asymptotically for large energies-allow us to describe the experimental data by the KLOE-2 collaboration [25], fulfill constraints from chiral perturbation theory on M C 1 , and reduce the freedom in the C-violating amplitudes to only one single complex normalization constant each. The phase of the latter is fixed by hermiticity and T violation, resulting in one real-valued free parameter for the isoscalar and isotensor transition, respectively. Ensuring that the Standard-Model contribution is in good accordance with the dispersive representation of Ref. [47], we extracted the contributions of M C 0 and M C 2 , whose interference with M C 1 give rise to the breaking of mirror symmetry in the η → π + π − π 0 Dalitz-plot distribution. We confirmed that the currently most precise measurement of the latter restricts the C-violating effects to a relative per mille level. Due to the strong kinematic suppression of M C 0 -the corresponding operator is smaller by four orders in the chiral expansion compared to M C 2 -the accompanying effective coupling constant g 0 is far less rigorously constrained than g 2 , by about three orders of magnitude.
Although there is no sufficiently precise Dalitz-plot measurement for η ′ → π + π − π 0 yet, we have demonstrated that, in principle, the larger available phase space lifts the suppression of the isoscalar C-odd amplitude to a large extent, making a potential mirrorsymmetry breaking therein more sensitive to M C 0 by roughly two orders of magnitude than in η → π + π − π 0 . Both decays would most likely be driven by the same, fundamental, BSM operators.
In a similar manner, we established a framework to analyze the decay η ′ → ηπ + π − , which is sensitive to another class of C-and CP -violating operators with isospin I = 1. In this decay the amplitude decomposes into the isoscalar SM amplitude M C 0 and a C-violating isovector contribution M C 1 . A regression to the Dalitz plot of the BESIII collaboration [49] yields again no evidence for C-violating effects and limits their patterns to a relative per mille level.