ω→3π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \rightarrow 3\pi $$\end{document} and ωπ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \pi ^{0}$$\end{document} transition form factor revisited

In light of recent experimental results, we revisit the dispersive analysis of the ω→3π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \rightarrow 3\pi $$\end{document} decay amplitude and of the ωπ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \pi ^0$$\end{document} transition form factor. Within the framework of the Khuri–Treiman equations, we show that the ω→3π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \rightarrow 3\pi $$\end{document} Dalitz-plot parameters obtained with a once-subtracted amplitude are in agreement with the latest experimental determination by BESIII. Furthermore, we show that at low energies the ωπ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \pi ^0$$\end{document} transition form factor obtained from our determination of the ω→3π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \rightarrow 3\pi $$\end{document} amplitude is consistent with the data from MAMI and NA60 experiments.


Introduction
A precise description of the amplitudes involving three particles in the final state is one of the open challenges in hadron physics. It becomes even more important in view of the high precision data from the GlueX, CLAS12, COMPASS, BESIII, and LHCb experiments, where various exotic states decaying to three particles have been or will be measured [1][2][3][4][5][6][7]. A proper description of three-particle amplitudes is also required for extraction of resonance parameters from lattice QCD computations [8][9][10][11][12].
Because of Bose symmetry only odd angular momentum is allowed in each of the ππ channels, and thus the final state is dominated by the J = I = 1 isobars, i.e. the ρ meson. The latter is related to the J = I = 1 ππ partial wave amplitude, which is known to high precision from the Roy analysis of ππ scattering [28]. The existing analyses of the decays of ω and φ to three pions [21,22,[29][30][31] are mainly based on unsubtracted dispersion relations, which result in a parameter-free prediction of the shape of the Dalitz plot. The φ decay was favorably compared to the high-statistics Dalitzplot data from the KLOE [32] and CMD-2 [33] experiments. Until recently the only available data for ω → 3π came from the WASA-at-COSY experiment [34]. Given that the nominal ρπ threshold is above the mass of the ω, the distribution of events in the Dalitz plot is rather smooth and, therefore, it can be efficiently parametrized in the experimental analyses by a low-order polynomial in the Dalitz-plot variables. The coefficient of the leading term in the Dalitz-plot polynomial expansion (i.e. the Dalitz-plot parameter α) obtained this way is consistent with the dominance of the ρ peak, even though it lies outside the kinematical boundary. However, the experimental uncertainties in the WASA-at-COSY measurement were too large to verify the prediction of dispersion relation calculations.
The situation changed recently when the high-statistics data from BESIII became available [35]. A new set of ω → 3π Dalitz-plot parameters was extracted from the data and was found to differ significantly from the predictions based on (unsubtracted) dispersion relations of the KT equations [21,22]. This is particularly unsettling since, as mentioned above, there is good agreement between the data and the predicted shape of the Dalitz plot in the case of the φ decays [21]. Therefore, in this paper we reanalyze the ω → 3π BESIII data with the KT equations.
As in any low-energy effective theory, the contribution from inelastic channels enters as (free) parameters, and this can be the origin of the discrepancy between the data and the calculation based on the unsubtracted dispersion relations. We also reanalyze the ωπ 0 transition form factor (TFF), which controls the ω → π 0 γ ( * ) amplitude. At low energies, the TFF is sensitive to the ω → 3π amplitude. There are recent data on the form factor from the MAMI [36] and NA60 [37,38] collaborations. As we will see below, the simultaneous analysis of both reactions allows one to better constrain the subtraction constant of the ω → 3π amplitude.
The analysis presented in this paper could also be relevant to understand the hadronic contributions to the anomalous magnetic moment of the muon. The presently observed ∼ 3σ deviation between theory [39][40][41][42] and experiment [43] has a potential to become more significant once new measurements at both Fermilab [44,45] and J-PARC [46] become available. The theoretical uncertainties mainly originate from the hadronic vacuum polarization (HVP) and the hadronic light-by-light (HLbL) processes, with ω → 3π, π 0 γ * amplitudes contributing to both. In particular, it was found that the reaction γ * → 3π , which builds upon V → 3π [47][48][49], gives the second-largest individual contribution to the HVP integral [50]. The same reaction together with the electromagnetic pion form factor constrain the doubly virtual pion transition form factor at low virtualities [48,49], which in turn gives the leading contribution to the HLbL process [51]. Additionally, in the calculation of the helicity partial wave amplitudes for γ * γ * → ππ [52][53][54], which are responsible for the twopion contribution to HLbL, the most important left-hand cut beyond the pion pole is almost exclusively attributed to the ω exchange. Thus, it depends on the ω → π 0 γ * TFF. Given the importance of the ω → 3π, π 0 γ * amplitudes and the fact that the currently available ones appear to be at odds with the most recent data, we find it timely to perform a new study of these reactions.
The paper is organized as follows. In Sect. 2 we briefly review the KT formalism for the ω → 3π decay, and show its relation to the ωπ 0 TFF. In Sect. 3 we discuss fits to the BESIII, MAMI and NA60 data. Our conclusions are given in Sect. 4. Details of the statistical analysis performed to determine uncertainties of the fits are given in the Appendix A.

Kinematics and initial definitions
We start by introducing the kinematical definitions for the ω( p V ) → π 0 ( p 0 ) π + ( p + ) π − ( p − ) process. The Mandelstam variables are defined as: Throughout this manuscript we work in the isospin limit with m 2 π = m 2 π ± = m 2 π 0 . The scattering angle in the s-channel, defined by the center of mass of the π + π − pair, is denoted by θ s and it is given by: where the momenta p(s) and q(s), are those of the π ± and π 0 , respectively, in the s-channel. The well-known Källen or triangle function λ(a, b, c) is defined as [55]: The also well-known Kibble function φ(s, t, u) is given by [56]: 5) and it defines the boundaries of the physical regions of the process through the solutions of φ(s, t, u) = 0. The Dalitzplot boundaries in t for a given value of s for the ω → 3π decay process lie within the interval [t − (s), t + (s)], with while the allowed range for s is:  (2.8) and at the same time decomposed into partial wave amplitudes where μναβ is the Levi-Civita tensor, μ ( p V , λ) is the polarization vector of the ω meson with helicity λ, and d J λ0 (θ s ) are the Wigner d-functions with θ s given by Eq. (2.2). For a V → 3π decay, H 0 = 0 and H + = H − , due to parity. As discussed in Ref. [22], one can rewrite the partial wave expansion for the invariant amplitude F(s, t, u) in the following form (2.10) where the exact relation between h (J ) + (s) and the kinematicsingularity-free isobar amplitudes f J (s) can be found in Ref. [22]. The KT representation of the invariant amplitude F(s, t, u) in Eq. (2.10) consists in substituting the infinite sum of partial waves in the s-channel by three finite sums of so-called isobar amplitudes, one for each of the s-, t-and u-channels. By truncating each sum at J max = 1 we obtain the crossing-symmetric isobar decomposition [21,22,26]: where each isobar amplitude, F(x), has only right-hand or unitary cut in its respective Mandelstam variable. For the ππ scattering a similar decomposition is known as the reconstruction theorem [57][58][59]. For J = 1, the relation between F(s) and f 1 (s) is obtained by projecting Eq. (2.11) onto the s-channel partial wave, where z s = cos θ s , and where the so-called inhomogeneitŷ F(s) contains the s-channel projection of the left-hand cut contributions due to the t-and u-channels. Its evaluation in the decay region requires a proper analytical continuation [60]. Assuming elastic unitarity with only two-pion intermediate states, we arrive at the KT equation for the ω → 3π decay, i.e. the unitarity relation for the isobar amplitude F(s): where δ(s) is the P-wave ππ phase shift. Given the discontinuity relation Eq. (2.14), one can write an unsubtracted dispersion relation (DR) for F(s) as which can be solved numerically [14,[19][20][21]. Its solution is given in terms of the usual Omnès function [61], defined by the real phase shift δ(s). For the latter we take the solution of the Roy equations of Ref. [28], that are valid roughly up to 1.3 GeV. From 1.3 GeV on we smoothly guide δ(s) to π to obtain the expected asymptotic 1/s fall-off behavior for the pion vector form factor (see e.g. Ref. [62]). The solution of Eq. (2.15) is written as: where the (complex) normalization constant a = |a| e iφ a is an overall normalization of the amplitude and can be factored out. Using PDG data, |a| can be fixed to reproduce the experimental ω → 3π decay width. No observables of the decay are sensitive to the overall phase φ a . Due to the asymptotic behavior of F(s) implied by Eq. (2.17), the amplitude F(s, t, u) satisfies the Froissart-Martin bound [21,63,64]. We emphasize that, even though F(s)/ (s) in Eq. (2.17) looks like a once-subtracted dispersion relation, F(s) actually satisfies the unsubtracted dispersion relation given in Eq. (2.15). Therefore, the energy dependence of F(s) is a pure prediction, which in the elastic approximation is given solely by the P-wave ππ phase shift. Note that Eq. (2.17) can be written in the form , (2.18) if b satisfies the following sum rule [21]: In Ref. [21] its value was computed, with the result: which we reproduce as a numerical cross-check. We note that, due to the three-particle cut, which become physically accessible in the decay amplitude, this subtraction constant is complex and is thus determined by two parameters, its modulus and phase, b = |b| e iφ b . In contrast to the unsubtracted DR in Eq. (2.15), one can start from a once-subtracted DR: The solution to Eq. (2.21) can be constructed as the linear combination [19,21]: where now b is not constrained to satisfy Eq. (2.19), and the functions F a (s) and F b (s) are given by These functions only need to be calculated once, since they are independent of the numerical values of a and b, which become fit parameters, as will be discussed in Sect. 3. For completeness, in Fig. 1 we show the solutions for F a (s) and F b (s) using a numerical iterative procedure similar to those employed in previous works [19,22,65,66]. By introducing one subtraction we reduce the sensitivity to the unknown high energy behavior of the phase shift and/or to the inelastic contributions, which are thus embedded in the subtraction constant. Furthermore, the parameter b allows to parametrize some unknown energy dependence of the ω → 3π interaction not directly related to ππ rescat-tering. 1 Strictly speaking, the amplitude F(s, t, u) built from F(s) in Eq. (2.22a) would not satisfy the Froissart-Martin bound [21,63,64] for an arbitrary value of the parameter b = b sum [cf. Eq. (2.19)]. In practice, however, given the lowenergy regime in which Eq. (2.22a) is applied, this bound is not relevant and we therefore do not constrain the value of b.
Finally, the measured differential decay width can be written in terms of the invariant amplitude F(s, t, u) as The ω → 3π Dalitz plot distribution is conventionally parametrized in terms of the variables X , Y defined by where s c = 1 3 (m 2 ω + 3m 2 π ) and R ω = 2 3 m ω (m ω − 3m π ). The {X, Y } variables are related to the polar ones {Z , ϕ} through X = √ Z cos ϕ and Y = √ Z sin ϕ, which enter into the Dalitz-plot expansion as: In Eq. (2.25), α, β and γ are the real-valued Dalitz-plot parameters and N is an overall normalization. In order to obtain α, β, and γ for a given theoretical amplitude F th (z, φ) we minimize [21] where D is the area of the Dalitz plot, φ(Z , ϕ) is φ(s, t, u) with s, t, and u expressed in terms of the polar variables {Z , ϕ}, and ξ 2 Dalitz denotes the average deviation of the theoretical description and the polynomial one relative to the Dalitz plot center. We also note that the Dalitz-plot parameters enter into the difference in Eq. (2.26) linearly, and thus the minimization can be algebraically solved.

General approach
The two amplitudes defined in the previous section depend on a total of five real parameters. The ω → 3π amplitude depends on |a| and b = |b| whereas the ωπ 0 transition form factor additionally depends on the subtraction constant at s = 0, f ωπ 0 (0), also complex.
For each of these sets we define the following χ 2 functions, where in χ 2 A2 and χ 2 NA60 the sum runs over the experimental points with √ s i 0.65 GeV. To determine the role of each data set, we start by considering the ω → 3π Dalitz plot parameters alone, since they only depend on |b| and φ b . In a first step, we fix |b| and φ b from the Dalitz plot parameters (i.e., by minimizing χ 2 DP ), and, in a second step, we fix |a| and f ωπ 0 (0) from the decay widths (i.e., by minimizing χ 2 ). We obtain the following values for the "2 par." case: whereas, for the "3 par." case, one gets: Because we are fitting two or three experimental points with two free parameters, the χ 2 DP is zero for the "2 par." case, and almost zero for the "3 par." case. In turn, this manifests in the large value of the errors shown in Eqs. (3.2). These errors are obtained through the condition χ 2 DP 1. We note that the value obtained for b is quite different from the value of b sum [cf. Eq. (2.20)], as also shown in Fig. 3. This reinforces the idea that, in order to achieve a proper description of the BESIII Dalitz plot parameters, an additional subtraction is needed within the KT formalism. 2 In Eqs. (3.2) we have fixed all the free parameters but φ ωπ 0 (0), and we now study the dependence of the ωπ 0 TFF on this phase. In Fig. 4 we show how χ 2 A2 and χ 2 NA60 depend on this phase for fixed values of the other parameters. We present the result for the "3 par." case, Eq. (3.2b), but an analogous result is obtained for the "2 par." case. We observe that there are two minima, one at φ ωπ 0 (0) 0.2 and another one at φ ωπ 0 (0) 2.5, to which we refer in what follows as "low φ ωπ 0 (0)" and "high φ ωπ 0 (0)" solutions, respectively. Furthermore, it is observed that the values of χ 2 NA60,A2 are similar in both cases, i.e., both solutions describe the data with similar quality.

Global fit results
Given that we are able to separately reproduce the experimental data on the two reactions, in the next step we perform a simultaneous fit. To that end, we minimize the following χ 2 -like function,   Table 2, case "3 par.". We show the results for the "low φ ωπ 0 (0)" (blue solid) and "high φ ωπ 0 (0)" (red dash-dotted) solutions. The error bands are obtained from a MC analysis of the fitted data, and represent the correlated 1σ uncertainty in our parameters where N DP = 2 or 3 is the number of Dalitz plot parameters considered, N = 2 the experimental partial widths, N A2 = 14 and N NA60 = 22 the experimental points in the two sets for |F ωπ (s)| 2 , and N = N DP + N + N A2 + N NA60 . This ensures that χ 2 functions with a smaller number of points are well represented in χ 2 , and are not overridden by those with a larger number of points. When the simultaneous fit is performed we observe, as expected, that the two solutions remain. The two minima are well separated, as can be seen in Fig. 4, so that we can analyze each solution individually. Besides these two solutions, we must also consider the two different sets of Dalitz-plot parameters given by the BESIII collaboration, as shown in Table 1. Therefore, we perform four different fits, and the fitted parameters, as well as the individual values of the χ 2 functions, are compiled in Table 2. We quote two errors, the first is statistical and the second one systematic/theoretical. The statistical errors are obtained through a Monte Carlo (MC) analysis with data resampling (bootstrap [72][73][74]), and they represent 1σ level uncertainties (see Appendix A for further details). The values obtained for the individual χ 2 functions imply a good quality of the fits. As a consistency check between the "2 par." and "3 par." data sets, we note that the values of the parameters are similar among the two "low φ ωπ 0 (0)" solutions (second and fourth columns in Table 2), as well as among the two "high φ ωπ 0 (0)" solutions (third and fifth columns). As an illustration, we show in Fig. 5 the function F(s) obtained using the values of the parameters that correspond to the "3 par." set, for both solutions. Regarding specifically the values of |b| and φ b , we note that both solutions fall well within the region determined by the fit to only BESIII data described in Sect. 3.1, see Fig. 3. This means that both solutions originate from that, but have much more constrained uncertainties as a result of the inclusion of the TFF data. We also note that the two widths considered in the χ 2 ( ω→3π and ω→π 0 γ ) are reproduced with the same central values and errors as the experimental ones.

3)
The systematic uncertainty attached to our calculation is obtained as follows. Because of the additional subtraction performed in the ω → 3π KT amplitude [Eqs. (2.22)] and the subtraction in the TFF [Eq. (2.28)], the convergence of the integral term in the latter is only as good as that of the case in which both amplitudes are obtained from unsubtracted dispersion relations. This convergence is relatively slow, so, Fig. 6 Dependence of the Dalitz plot parameters (X = α, β, γ ) when a cutoff is used as an upper integration limit in the TFF dispersive representation, Eq. (2.28). This dependence is used to estimate a systematic uncertainty. See the text for more details in order to quantify potential uncertainties from it, we study the dependence of our fits for a finite upper limit in the integral in Eq. (2.28). In particular, we fit the free parameters for running finite values of the cutoff , obtaining also the output quantities (Dalitz plot parameters and TFF) for these cutoff values. In general, we observe that the results for a cutoff 1 GeV 2 are already very similar to the case when = ∞. For instance, in Fig. 6 we represent the dependence on of the Dalitz plot parameters (to be discussed later on), for the "3 par.", "low φ ωπ 0 (0)" case. The uncertainty we attach to every quantity is the absolute value of the difference between the fits performed for = ∞ and = (m ω + m π ) 2 . This latter value is the ωπ threshold, beyond which our model would be certainly not reliable. We observe that the uncertainty attached to the fitted parameters and to the Dalitz plot parameters is sizeable, but not for the TFF itself, due to the fact that it dominates the χ 2 -like function. For this reason, we do not show a systematic uncertainty of the TFF in what follows.
The results for the TFF are shown in Fig. 7 for the low and high φ ωπ 0 (0) solutions. It can be seen that both of them agree very well with the experimental points, except for the highest two points of the NA60 data. 3 Also, it should be noted that both solutions are almost indistinguishable. The largest difference is at the ωπ 0 invariant mass √ s 0.3 GeV, which is near the 2π threshold, but even there they are compatible at 1σ level. Although we will later on compare in detail our results with other approaches, it is worth pointing out 3 These two points give a contribution of around 17 to χ 2 NA60 . However, we note that fits without these two points give similar results as the ones discussed in the text.

Fig. 7
Normalized ωπ 0 TFF, | f ωπ 0 (s)/ f ωπ 0 (0)| 2 . The data are taken from Refs. [36][37][38]. The lines, and their associated error bands, represent our two different solutions, which overlap almost completely in the ωπ 0 invariant mass range shown. The case shown here is that corresponding to the "3 par." fit. The curves for the "2 par." case are very similar. For comparison, we also show the Vector Meson Dominance (dot-dotdashed brown) prediction, and that of the model without KT equations (dotted pink curve) discussed in Sect. 3.3, cf. Eqs. (3.4) here that our theoretical description of the data represents an improvement over previous theoretical analyses [22,68,75].
We note that the different phase φ ωπ 0 (0) in both solutions translates into a difference in the phase of the TFF in a large region of ωπ 0 invariant mass, up to √ s 0.6 GeV, as shown in Fig. 8. For energies √ s 0.6 GeV the phase motion associated with the ρ meson kicks in, and both solutions approximately converge. This phase, or more properly the phase difference φ ωπ 0 (0) − φ a (see Sect. 2.3) has not been measured, to the best of our knowledge, and thus Fig. 8 constitutes a prediction for it. 4 Finally, we note that the "low φ ωπ 0 (0)" solution is rather close to φ ωπ 0 (0) = 0, and the "high φ ωπ 0 (0)" is close (but less than the previous one) to φ ωπ 0 (0) = π . If the amplitudes were computed from a Lagrangian approach with a stable ω, the couplings in the Lagrangians would be real. Then, one would expect real values for a and f ωπ 0 (0), and thus their relative phase could only be 0 or π . Anyhow, we find that the inclusion of this phase with a value different from 0 or π improves the description of the data, since they are different from zero by approximately 2σ .
In what relates to the Dalitz plot parameters, we find good agreement between the input taken from BESIII and our results, see Table 1, which results in the low χ 2 DP shown in Table 2, in the four cases considered (low or high φ ωπ 0 (0), 2 or 3 Dalitz plot parameters). The largest difference between  Table 2. The error bands represent our (correlated) 1σ uncertainties in the fitted parameters, obtained from a MC analysis of the data. We show here the curves for the "3 par." fit. The phases for the "2 par." case are very similar. For comparison, we also show the prediction of the model without KT equations (dotted pink curve) discussed in Sect. 3.3, cf. Eqs. (3.4) observables used in our fit for the "3 par." case is found in γ . The values that we obtain, γ = (19 ± 5 ± 4) × 10 −3 and (29 ± 6 ± 8) × 10 −3 for the "low φ ωπ 0 (0)" and "high φ ωπ 0 (0)" solutions, respectively, are both compatible with the experimental one used in the fit, γ = (22 ± 29) × 10 −3 . However, our values are found to be better constrained and indicate that this parameter is non-zero at a ∼ 3σ level. Interestingly, the two values of γ are only marginally compatible and a more precise measurement of the ω → 3π Dalitz-plot parameters could help in pinning down the correct solution. A similar argument, though less stringent, can be made for β in the "2 par." fits. As previously mentioned, in the process of obtaining the systematic errors described above we have also recomputed the Dalitz plot parameters, which allows us to propagate this uncertainty into them. This uncertainty is also shown in Table 1, and turns out to be small for α and β, but larger for γ .

Comparison with previous approaches
Our results obtained by solving KT equations for the ω → 3π amplitude, are compared with those from Refs. [21,22] in Table 1. The difference between these approaches and ours lies in the subtraction that we have performed on the KT dispersion relations, which introduces an additional free param-eter, b. In Ref. [21], an estimation for this parameter is given by enforcing the once-subtracted DR to be equivalent to the unsubtracted DR. This value, b sum 0.55e 0.15i GeV −2 , Eq. (2.20), turns out to be far away from our fitted b (for any of the fits in Table 2), which reaffirms the need of the extra subtraction. Due to this subtraction, and the fits performed in Sects. 3.1 and 3.2, our results for the Dalitz-plot parameters are in agreement with those of the BESIII experiment. Nonetheless, we must recall here that the disagreement between the BESIII [35] Dalitz-plot parameters and that obtained with unsubtracted KT equations [21,22] lies only in the α parameter, as can be seen in Table 1. In particular, the 3-parameter determinations give α = 111(18), 77, and 80, respectively. Our results (with a once-subtracted KT equation) for the α parameter are 109 (14)(2) and 112 (15)(2) for the high and low-φ ωπ 0 (0) solutions, respectively. While one subtraction is required to achieve good agreement with the experimental data, the difference between the oncesubtracted and unsubtracted KT equations determination of the α parameter turns out to be only around 2σ . 5 The values of the TFF given by the KT approach without the additional subtraction used in our work for the ω → 3π amplitude lie systematically below the experimental points [22,68]. In Ref. [22] it was shown that without the extra subtraction a satisfactory result for the TFF can be obtained only if additional terms are retained in the non-dispersive term (see Fig. 8 of that reference). In contrast, as discussed in Sect. 3.2, our results for the TFF are in good agreement with the experimental data. In particular, our approach represents a significant improvement in the description of the NA60 higher energy data, although we will further discuss about this issue in Sect. 3.4.
In Table 1 we show the results obtained in Refs. [21,22] when the crossed channel effects, which are the essential outcome of the KT equations, are "turned off" from the isobar F(s). In practical terms, this is achieved by neglecting the contribution ofF(s) in Eq. (2.22), such that F(s) is simply an Omnès function times a constant, The reduced full amplitude would then read The proportionality constant, a instead of a, is chosen to reproduce the ω → 3π width, 10 −2 a = 2.818(18) GeV −3 , but it is a global constant and does not affect the values of the Dalitz plot parameters. Interestingly, as discussed in Sect. 1, the Dalitz plot parameters obtained in Refs. [21,22] in this simplified approach appear to be in better agreement with the recent experimental determination Fig. 9 Absolute value (green solid), real (red dot-dot-dashed) and imaginary (blue dashed) parts of the F(s)/F sim (s) ratio, as described in the text. The ratio is shown in the physical decay range of s by BESIII [35] than those obtained with the crossed channel effects included (but no extra subtraction), cf. Table 1, rows denoted "w/o KT" vs. those denoted "w KT", respectively. In sharp contrast, we show in this work that the results we obtain by keeping the crossed channel effects, and with the additional subtraction, reproduce very well the experimental Dalitz-plot parameters, and are consistent with the ωπ 0 TFF. We first discuss why the determination of the Dalitzplot parameters is very similar in our approach (subtracted KT) and in the simpler model (no KT, Eqs. (3.4)). Later on, we will compare the results for the TFF.
The aforementioned agreement is clear, as can be seen in Table 1, and hence there must be some sort of cancellation that "brings back" our full subtracted KT approach into the simpler, no KT model. Naively, if one thinks that the KT formalism is overestimating the crossed channel effects, it would be expected that this cancellation would occur in the isobar amplitude itself, F(s), i.e., that the effect of the crossed channels is mostly linear and thus can be absorbed by the additional subtraction constant, b. In this case, the ratio F(s)/F sim (s) should be essentially constant. We show in Fig. 9 that this is certainly not the case, although the modulus of the ratio is still around 1. Here, we are taking the parameters of the "low φ ωπ 0 (0)" solution for the "3 par." case, but similar results are obtained in the other fits. This demonstrates that the cancellation is not trivial, as one would expect if the crossed channel effects were simply being overestimated.
The cancellation must thus occur at the level of the squared amplitude, |F(s, t, u)| 2 = |F(s) + F(t) + F(u)| 2 . In Fig. 10 we show |F(s, t, u)| 2 for √ s in the physical decay region for two lines across the (s, t, u) plane, namely, t = u and t = s c (respectively corresponding to X = 0 and Y = − √ 3X in the usual (X, Y ) Dalitz plot variables, cf. Eq. (2.24)). We also show in the figures the function |F sim (s, t, u)| 2 , i.e., the full amplitude squared for the simpler model [Eqs. (3.4)] discussed above. It can be seen in Fig. 10 that the differences between both squared moduli are quite small. We also find that the phase difference between our F(s, t, u) and F sim (s, t, u) is essentially constant. This large cancellation explains the coincidence of the results for the Dalitz plot parameters in both approaches. 6 As a result of the above discussion, one might question the necessity of the full approach if, after all, the rather simpler description with no subtractions and no crossed channel effects, Eq. (3.4), seems to work just fine. However, it must be noted that this simpler model only describes well the ω → 3π Dalitz-plot parameters, but not the distribution for the φ → 3π decay [21] nor the more precise experimental information on the ωπ 0 TFF. In Ref. [68] it is shown that a model which ignores the crossed-channel effects by inserting f 1 (s) = a (s) into Eq. (2.28) gives a result well below the experimental points (see Fig. 5 of Ref. [68]). We could also take the partial wave that results from Eqs. (3.4a) and (3.4b), which is given by This model, when introduced into Eq. (2.28), produces the result shown as a pink dotted line in Fig. 7, which is well below the experimental points and our results. 7 This result for the TFF is very similar to that of Ref. [68] mentioned above.
In summary, from a phenomenological point of view, our description of the Dalitz-plot parameters and of the TFF using a once-subtracted version of the KT equations is (not surprisingly) better than that obtained with unsubtracted KT equations [21,22,68]. On the other hand, the simpler model of Eqs. (3.4), in which the KT effects are ignored, describe properly the Dalitz-plot parameters (see the discussion above about Figs. 9 and 10), but not the TFF data. Therefore, it seems that our approach, in which a KT equation for the ω → 3π amplitude is solved with an additional subtraction, is the minimal theoretical setup that is able to simultaneously describe both sets of data. From a more theoretical perspective, it is clear that the crossed channel effects must be present in any 2 → 2 or 1 → 3 amplitude, even if they are negligible or can be mimicked by polynomial terms [59]. The KT formalism offers a simple framework which allows to provide the partial waves in the direct channel with left hand cuts in terms of the isobars of the crossed channels, while In Refs. [29,30] the authors derived bounds on the modulus of the TFF based on unitarity and analyticity using the methods of unitary bounds [78,79]. The upper and lower bounds of the TFF are given by: where: 9 The conformal mapping z(s), maps the elastic region 4m 2 π s s + = (m ω + m π ) 2 into the real segment z(4m 2 π ) z(s) 1, whereas the inelastic region s s + is mapped into the disk |z(s)| = 1. The so called outer function C(z) is given by: The quantities I and g(0) are written as: (3.9c) As can be seen from the previous discussion, the central value of the bounds f (0) ωπ 0 (s) depends only on input quantities, i.e., the discontinuity disc f ωπ 0 (s) for s s + and the subtraction constant f ωπ 0 (0). In Fig. 11 we show as a green, solid line this central value for the case in which the unsubtracted ω → 3π KT equations [22,68] are used as an input for the discontinuity. For comparison, we also show the TFF itself with a green, dashed line, and it can be seen that both curves are similar. In contrast, the amplitude of the bound δf (0) ωπ 0 (s) is proportional to I , Eq. (3.9b), which involves the constant I , Eq. (3.9c). This constant, in turn, depends on the form factor at high energies, s > s + . In Ref. [30] the value I = 4.63 GeV −4 is estimated by interpolating the high energy data by the CMD Collaboration [80], which come from the e + e − → ωπ 0 reaction. They are shown in Fig. 11, together with those of the SND Collaboration [81] (see also Refs. [82,83]). These data are measured at energies above the ωπ 0 threshold, which lies beyond the focus of our paper. This value I = 4.63 GeV −4 is the one we have used to compute δf (0) ωπ 0 (s), and in Fig. 11 we show with a green band the resulting range f ± ωπ 0 (s) for the unsubtracted ω → 3π amplitude. Now we discuss the same calculation of the bounds when the once-subtracted ω → 3π amplitude is used. In particular, we use the 3 par., low φ ωπ 0 (0) fit in Table 2 (similar  results are Fig. 11 with a blue, solid line. The blue, dashed line shows the TFF (i.e., the blue line of Fig. 7), and it can be noted that, as in the case of the unsubtracted ω → 3π amplitude, both curves are similar. However, we cannot compute the amplitude of the bound δf (0) ωπ 0 (s) because if we use the previous value of I into our calculation of I , the r.h.s. in Eq. (3.9b) turns out to be negative. This is due to the fact that the integral terms in Eq. (3.9b) are much larger when the once-subtracted ω → 3π amplitude is used instead of the unsubtracted one. In Ref. [30] the possibility of Eq. (3.9b) being negative is interpreted as an incompatibility of the computed form factor with the principles of unitarity and analyticity. From our perspective, since the value I = 4.63 GeV −4 is obtained from higher energy experimental data, the incompatibility is actually between these data and our low energy description of the transition form factor. Since, ultimately, our low energy TFF is driven by the low energy data, this incompatibility could actually point to a tension between low and high energy data. No strong claims can be done in this respect since an extension to the higher energy region would require taking into account higher resonances in the pion vector form factor, and possibly the inclusion of inelastic thresholds. Still, in order to give an estimation of the bounds in Eq. (3.5), instead of taking the value I = 4.63 GeV −4 from Ref. [30] we can compute I in Eq. (3.9c) from our own model (we use the 3 par., low φ ωπ 0 (0) solution of Table 2). In this way, we obtain a larger value for I , that we denote I th , I th = 33 GeV −4 , and now the quantity I 2 in Eq. (3.9b) becomes positive, "confirming" the compatibility of our form factor with analyticity and unitarity. If we use this value I th in the calculation of δf (0) ωπ 0 (s) for the once subtracted ω → 3π amplitude, we obtain the blue band shown in Fig. 11. The band is larger than in the unsubtracted case because of the larger value of I , as discussed in Ref. [30]. Interestingly enough, we see that even if the band is larger, the two highest energy data points of the NA60 collaboration still lie outside of the bounds, pointing to the aforementioned possible tension. We thus conclude that more analysis is needed, both on the theoretical and on the experimental sides, to bring together the low and high energy data and their theoretical description.

Outlook
Summary In this work we have explored the benefits of a simultaneous analysis of the ω → 3π decay and the ωπ 0 transition form factor. The motivation for this study is manifold. First, from the point of view of strong interactions, the decay ω → 3π offers a good environment to study the dynamics of the ππ subsystems under rather clean conditions. Second, the BESIII collaboration has reported a highstatistics measurement of the ω → 3π Dalitz plot distribution, and pointed out a possible overestimation of the crossedchannel contributions in the KT equations. Third, there are recent data on the shape of the ωπ 0 TFF from the MAMI and NA60 collaborations making such an analysis of timely interest.
For the ω → 3π amplitude we follow a dispersive representation with subtractions that emerges from the solution of the KT equation. It thus satisfies the constraints posed by analyticity (to some extent), crossing symmetry and (elastic) unitarity, and it is completely determined by the ππ P-wave scattering phase shift, except for the values of the subtraction constants. In this work we have performed one subtraction, which introduces an additional free parameter, b, apart from the usual global normalization a that is fixed from the partial decay width. We fix this extra parameter, which is characterized by its modulus |b| and phase φ b , from fits to experimental data. The ω → 3π amplitude, in turn, enters the once-subtracted dispersive parametrization of the ωπ 0 TFF Eq. (2.28), introducing its phase at s = 0, φ ωπ 0 (0), as a new ingredient of this work.
Our first analysis proceeds in two steps. On a first step, we use the two different sets of Dalitz-plot parameters given by BESIII and the corresponding partial decay widths to fix all free parameters (|a|, |b|, φ b , | f ωπ 0 (0)|) except for φ ωπ 0 (0). These results bring us to a first relevant observation: the value of the subtraction constant b needed to faithfully reproduce the Dalitz-plot parameters is found to be significantly different (see Fig. 3) from the sum-rule value estimated from the unsubtracted version of the KT equations. On a second step, the dependence of the ωπ 0 TFF on φ ωπ 0 (0) is studied in relation to the MAMI and NA60 data. It is found that there are two well separated minima in this variable.
We have also performed a combined analysis to all available experimental information including Dalitz-plot parameters and form-factor data, and observed that the two solutions for φ ωπ 0 (0) remain. Interestingly enough, the values for the subtraction constant b obtained from the joint fits have a much better constrained uncertainty than that in the individual fits to the BESIII Dalitz-plot parameters (see Fig. 3), however being in perfect agreement with it. This reaffirms the need of the additional subtraction constant, although the difference in the determination of the α Dalitz-plot parameter in the unsubtracted [21,22] and once-subtracted KT equations is around 2σ .
From the Dalitz-plot parameters associated to our combined fits (see Table 1), we can draw a second relevant observation. While the values that we obtain for the Dalitz-plot parameters are found to be in agreement with the experimental ones, our values carry a smaller error and indicate a significance for the Dalitz-plot parameter γ of ∼ 3σ , including statistical and systematic errors. The former turn out to be the dominating ones in α and β, while the latter is more relevant for γ . Furthermore, our results for the normalized ωπ 0 TFF (Fig. 7) show a satisfactory description of the experimental data, except for the highest two points of the NA60 collaboration.
Open questions Even though we achieved a simultaneous description of the Dalitz-plot parameters and the TFF data, it comes as a surprise that the predictions for the ω → 3π amplitude are so different between the unsubtracted and once-subtracted versions of the KT equations. (This can be visualized either in the discrepancy between the Dalitz-plot parameters in both cases, or in the large difference between the fitted subtraction constant b respect to the sum-rule expectation.) Moreover, this does not seem to happen in φ → 3π , despite the larger phase space, which makes this difference even more intriguing.
It is also important to note that, due to the goal of our work, the analysis of the ωπ 0 TFF has been restricted to the relatively low energy region of the NA60 (ω → π 0 μ + μ − ) and MAMI (ω → π 0 e + e − ) data. Because of this, we have not taken into account the higher energy region beyond the ωπ 0 threshold, where there are experimental data [80][81][82][83] coming from the reactions e + e − → ωπ 0 . To do so would require to consider also higher resonances in the ππ phase shifts and the inclusion of inelastic thresholds, something clearly outside the scope of the present analysis. Furthermore, the NA60 data currently have much smaller uncertainties than the MAMI ones, which translates into the fact that our fits to the TFF have been dominated by the former, with almost no influence of the latter. The NA60 data drive the TFF curve towards higher values (even more if one aims to describe also the last two NA60 data points). While a fit of an additional subtraction constant is able to reduce the tension with the NA60 data, this certainly influences the extrapolation to higher energies. In particular, we have seen in Sect. 3.4 that this behaviour impacts the calculation of the bounds of the TFF discussed in Refs. [29,30]. In turn, this could point out to inconsistencies with the CMD and SND scattering data, a conclusion similar to those of Refs. [29,30]. We have also seen that the experimental points of the NA60 collaboration in the range √ s = 0.6−0.7 GeV lie outside of these bounds. Therefore, we hope that our study strengthens the case for a reanalysis of all these decays and/or new measurements thereof, either to reduce uncertainties or to address eventual incompatibilities.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study and no experimental data has been listed.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Fig. 12
One dimensional distributions of the free parameters in the Monte Carlo analysis performed, for the "3 par." case fit. A Gaussian distribution with the average and error quoted in Table 2 is superim-posed for each solution. The double arrows on the upper part of each histogram represent the 1σ uncertainty intervals Fig. 13 One dimensional distributions of our computed Dalitz plot parameters in the Monte Carlo analysis performed, for the "3 par." case (top row) and "2 par." case (bottom row). A Gaussian distribution with the average and error quoted in Table 1 is superimposed for each solution. The double arrows on the upper part of each histogram represent the 1σ uncertainty intervals. The green line represents the BESIII experimental determination, also shown in Table 1 of the four fits considered in Table 2, we generate O(10 4 ) sets of the data (resampling) described in Sect. 3.1, each single datum following a Gaussian distribution. For each of these sets, a fit is performed and each of the output quantities of our work (DP parameters, TFF, etc.) are computed for that fit. In this way, all possible known correlations are taken into account. The values obtained in this work quoted in Tables 1  and 2, as well as those represented in Figs. 5, 7, and 8 are the average value and the standard deviation of those quantities in all the fits generated.
In the histograms of Figs. 12 and 13 we show the probability distribution of the fitted free parameters and the resulting contours computed with the average and standard deviation, and correlation parameter of each pair. This provides a pictorial representation of the correlation and accounts for the non-Gaussianity of the distributions Dalitz plot parameters, respectively, obtained in our MC analysis for both the low and high φ ωπ 0 (0) solutions. We show the "3 par." case, but similar results are seen for the "2 par." case. In general, the parameters are seen to follow a Gaussian distribution, although some deviations are seen from this behaviour, specially for |b| and φ b . This non-Gaussianity is, of course, inherited from the χ 2 DP function, as can be seen in Fig. 3.
The correlation parameter between the fitted parameters and/or the computed quantities can be calculated in a standard way. However, the two-dimensional distributions are not always Gaussian, and we therefore prefer to show the two-dimensional projections of (a small sample of) our MC simulations in Fig. 14.