Dispersion relations for
Abstract
We present a dispersive analysis of the decay amplitude for that is based on the fundamental principles of analyticity and unitarity. In this framework, finalstate interactions are fully taken into account. Our dispersive representation relies only on input for the and scattering phase shifts. Isospin symmetry allows us to describe both the charged and neutral decay channel in terms of the same function. The dispersion relation contains subtraction constants that cannot be fixed by unitarity. We determine these parameters by a fit to Dalitzplot data from the VES and BESIII experiments. We study the prediction of a lowenergy theorem and compare the dispersive fit to variants of chiral perturbation theory.
∎ \DeclareUnicodeCharacter00B0°
e1email: isken@hiskp.unibonn.de \thankstexte2email: kubis@hiskp.unibonn.de \thankstexte3email: schneider@hiskp.unibonn.de \thankstexte4email: pstoffer@ucsd.edu
1 Introduction
The treatment of hadronic threebody decays using dispersion relations is a classic subject. Already in the 1960s, Khuri and Treiman developed a framework in the context of decays Khuri:1960zz (). One of its main virtues is the fact that the most important finalstate interactions among the three pions are fully taken into account, in contrast to perturbative, fieldtheorybased approaches: analyticity and unitarity are respected exactly. This becomes the more important, the higher the mass of the decaying particle, hence the higher the possible energies of the twopion subsystems within the Dalitz plot. But even in decays of relatively light pseudoscalar mesons like , finalstate interactions strongly perturb the spectrum. In such a case, a dispersive approach that resums finalstate rescattering effects is essential to reach high precision; see Refs. Kambor:1995yc (); Anisovich:1996tx (); Kampf:2011wr (); Lanz:2013ku (); Guo:2015zqa (); Guo:2016wsi (); Colangelo:2016jmc (); Albaladejo:2017hhj (). In this article, we present the application of these techniques to the decay .
The decay has received considerable interest in past years for several reasons. Due to the U(1) anomaly the is not a Goldstone boson and therefore “standard” chiral perturbation theory (PT) based on the spontaneous breaking of SU(3)SU(3) chiral symmetry fails to adequately describe processes involving the . In the limit of the number of colors becoming large (“large limit”) the axial anomaly vanishes, which leads to a U(3)U(3) symmetry, so that a simultaneous expansion in small momenta, small quark masses, and large gives rise to a power counting scheme that in principle allows one to describe interactions of the pseudoscalar nonet (). However, the question whether this framework dubbed large PT Kaiser:2000gs (); Escribano:2010wt () is actually wellestablished remains under discussion, mainly due to the large mass. This is an issue that can in principle be addressed by a study of . So far there are indications that a large PT treatment alone is not sufficient to describe the decay, as finalstate interactions play a rather important role, see Refs. Escribano:2010wt (); Riazuddin:1971ie ().
Furthermore, the decay channel could be used to constrain scattering: the mass is sufficiently small so that the channel is not polluted by nonvirtual intermediate states other than the rather wellconstrained scattering. In the past claims were made that the mechanism via the intermediate scalar resonance even dominates the decay Deshpande:1978iv (); Singh:1975aq (); Fariborz:1999gr (). These claims are based on effective Lagrangian models with the explicit inclusion of a scalar nonet incorporating the , , and [] resonances. They were further supported by Refs. Beisert:2002ad (); Beisert:2003zs (): a chiral unitary approach shows large corrections in the channel and there is a dominant lowenergy constant in the U(3) PT calculation that is saturated mostly by the . The wave, however, was found to be strongly suppressed Borasoy:2005du (); RobinDiss (); Kubis:2009sb ().
The decay channel is expected to show a cusp effect at the chargedpion threshold Kubis:2009sb () that in principle can be used to obtain information on scattering lengths. So far this phenomenon has not been observed: the most recent measurement with the GAMS spectrometer did not have sufficient statistics to resolve this subtle effect Blik:2009zz ().
The extraction of scattering parameters such as the scattering length and the effectiverange parameter is a more complicated subject compared to scattering. There is no oneloop cusp effect as in the channel, since the threshold sits on the border of the physical region and not inside. The hope of extracting scattering parameters from a twoloop cusp is shattered likewise: there is a rather subtle cancellation of this effect at threshold (see Refs. Schneider:2009rz (); Diplom () for an elaborate discussion).
Measurements of the Dalitz plot of the charged channel have been performed by the VES Dorofeev:2006fb () and BESIII Ablikim:2010kp () collaborations, while earlier measurements at rather low statistics have been reported in Refs. Kalbfleisch:1974ku (); Briere:1999bp (). The more recent measurements seem to disagree considerably with regard to the values of the Dalitzplot parameters, and also in comparison with the GAMS measurement Blik:2009zz () the picture remains inconsistent.
This article is structured as follows. We will start by discussing the necessary kinematics as well as the resulting analytic structure of in Sect. 2, before deriving and analyzing dispersion relations for the decay in Sect. 3. In Sect. 4, we will discuss the numerical solution of the dispersion relation. The results of the fits to data will be discussed in Sect. 5. Predictions for higher Dalitzplot parameters, the occurrence of Adler zeros close to the softpion points, and predictions for the decay into the neutral final state are discussed in Sect. 6. Finally, we perform a matching of the free parameters to extensions of PT in Sect. 7. Some technical details are relegated to the appendices.
2 Kinematics
We define transition amplitude and kinematic variables of the decay in the usual fashion,
(1) 
where refer to the pion isospin indices.^{1}^{1}1In the following, we will consider both the charged decay channel and the neutral channel . They differ only by isospinbreaking effects. We define the Mandelstam variables for the threeparticle decay processes according to
(2) 
which fulfill the relation
(3) 
The process is invariant under exchange of the pions, that is, under . In the centerofmass system of the two pions, one has
(4) 
where refers to the scattering angle,
(5) 
with the Källén function and . Similarly, in the centerofmass system of the channel, one finds
(6) 
with and
(7) 
Due to crossing symmetry, the channel relations follow from , .
The physical thresholds in the three channels are given by
(8) 
3 Dispersion relations for
In this section, we set up dispersion relations for the decay process , in analogy to previous work on different decays into three pions Khuri:1960zz (); Anisovich:1996tx (); Kambor:1995yc (); Niecknig:2012sj (). The idea is to derive a set of integral equations for the scattering processes and with hypothetical mass assignments that make these (quasi)elastic: in such a kinematic regime the derivation is straightforward. The dispersion relation for the decay channel is then obtained by analytic continuation of the scattering processes to the decay region.
We will begin our discussion by decomposing the amplitude in terms of functions of one Mandelstam variable only. This form will prove very convenient in the derivation of the integral equations and their numerical solution at a later stage. Such a decomposition goes under the name of “reconstruction theorem” and was proven to hold in the context of chiral perturbation theory up to twoloop order for pion–pion scattering Stern:1993rg (). It was subsequently generalized to the case of unequal masses Ananthanarayan:2000cp () and to general scattering of pseudoscalar octet mesons Zdrahal:2008bd (). We derive it in A, finding the form
(9) 
where are functions of one variable that only possess a righthand cut. Here, refers to angular momentum and to isospin: isospin conservation of the decay constrains the total isospin of the finalstate pion pair to , while the system always has . Equation (3) follows from a partialwave expansion of the discontinuities in fixed, , and  dispersion relations, symmetrized with respect to the three channels. Given the smallness of the available phase space, the partialwave expansion is truncated after  and waves. A wave contribution is forbidden by charge conjugation symmetry. We stress that the truncation only neglects the discontinuities or rescattering phases in partial waves of angular momentum : projecting of Eq. (3) on the wave (in the channel) yields a nonvanishing result, however, this wave is bound to be real apart from threeparticlecut contributions.
We will briefly discuss the finalstate scattering amplitudes that are involved in , namely and . Given again the maximum energies accessible in the decay, both rescattering channels are treated in the elastic approximation, such that the corresponding partial waves can be parametrized in terms of a phase shift only, without any inelasticity effects. The scattering amplitude (confined to ) is approximated by
(10) 
with denoting the wave phase shift. Analogously, the scattering amplitude can be represented, neglecting  and higher waves, according to
(11) 
where is the phase shift of angular momentum .
The unitarity condition for the decay of the to a generic threebody final state can be written as
(12) 
where denotes the decay amplitude and describes the transition, while the sum runs over all possible intermediate states .^{2}^{2}2Here and in the following, relations that involve the discontinuity are always thought to contain an implicit function that denotes the opening of the respective threshold, i.e. for the channel and for the channel. The integration over the intermediatestate momenta is implied in this shorthand notation. Limiting the sum to and rescattering, carrying out the phasespace integration, and inserting the partialwave expansion for the and amplitudes Eqs. (10) and (3), we find
(13) 
where denotes the angular integration between the initial and intermediate state of the respective , , channel subsystem, while refers to the centerofmass scattering angles between the intermediate and final state. Finally, we can insert the decomposition of the decay amplitude (3) on the left and righthand side of Eq. (3) and find the unitarity relations for the singlevariable functions :
(14) 
where . The inhomogeneities are given by
(15) 
where we have defined the shorthand notation
(16) 
Note that the analytic continuation of Eqs. (3) and (3) both in the Mandelstam variables and the decay mass involves several subtleties. This is discussed in detail in Refs. Kambor:1995yc (); Anisovich:1996tx (); Walker (); StefanDiss () for , as well as in Ref. Niecknig:2012sj () for and in Ref. Schneider:2012nng () specifically for . One important consequence is the generation of threeparticlecut contributions in the decay kinematics considered here.
The solutions of the unitarity relations Eq. (3) can be written as
(17) 
where the Omnès functions are given as
(18) 
(with the appropriate thresholds or ). The order of the subtraction polynomials in the dispersion relations is determined such that the dispersion integrals are convergent. However, we can always “oversubtract” a dispersion integral at the expense of having to fix the additional subtraction constants and possible ramifications for the highenergy behavior of our amplitude.^{3}^{3}3Each additional subtraction constant contributes an additional power of to the asymptotic behavior of the amplitude if the corresponding sum rule for the subtraction constant is not fulfilled exactly. This can lead to a violation of the Froissart–Martin bound. To study the convergence behavior of the integrand we have to make assumptions as regards the asymptotic behavior of the phase shifts. We assume
(19) 
as . Note that an asymptotic behavior of implies that the corresponding Omnès function behaves like in the same limit.
Finally, we assume an asymptotic behavior of the amplitude inspired by the Froissart–Martin bound Froissart:1961ux (); Martin:1962rt (),
(20) 
which allows the following choice for the subtraction polynomials:
(21) 
The subtraction constants thus defined are correlated since the decomposition Eq. (3) is not unique. By virtue of Eq. (3), there exists a fourparameter polynomial transformation of the singlevariable functions that leaves invariant. Restricting the asymptotic behavior to Eq. (20) reduces it to the following twoparameter transformation:
(22) 
which can be used to set the first two coefficients in the Taylor expansion of around to zero. Since the transformation polynomial is a trivial solution of the dispersion relation (with vanishing discontinuity), the transformed representation still can be cast in the form of Eq. (3). Relabeling the transformed subtraction constants and inhomogeneities to the original names, we obtain the following form of the integral equations:
(23) 
In the following, we will neglect the wave as it is strongly suppressed with respect to the wave of and scattering; see for example Refs. Kubis:2009sb (); Bernard:1991xb (). In fact, the wave has exotic quantum numbers, such that the phase shift is expected to be very small at low energies. In chiral perturbation theory, this phase (or the corresponding discontinuity) only starts at (three loops) and is therefore, in this respect, as suppressed as all higher partial waves.
The decomposition of the amplitude in this case simply reads
(24) 
We will call the dispersive representation as outlined above “DR”, as it depends on four subtraction constants. In our numerical analysis, we compare it to a representation where we further reduce the number of free parameters by assuming a more restrictive asymptotic behavior of the amplitude: for large values of , , , respectively. In this case, the symmetrization procedure in the reconstruction theorem is possible for waves only and the singlevariable functions fulfill the integral equations
(25) 
Note that the transformation (3) still allows us to set the first subtraction constant in to zero. The second subtraction constant cannot be removed without changing the asymptotic behavior. As there are three subtraction constants , , and , we refer to this setup as “DR”.
The representation DR (3) (with ) is equivalent to DR (3), provided that the subtraction constants and fulfill a certain sum rule in order to guarantee the constraint of the asymptotic behavior. Explicitly, the relation between the DR and DR subtraction constants is given by
(26) 
where , , and the sum rule is encoded by the integrals
(27) 
Since the Froissart–Martin bound is strictly valid only for elastic scattering, a pragmatic approach concerning the number of subtractions is advisable. On the one hand, we try to work with the minimal number of subtractions allowing for a good fit of the data. On the other hand, additional subtractions help to reduce the dependence on the highenergy behavior of the phase shifts. Therefore, in our analysis we compare both representations DR and DR.
As we will show, the representation DR of Eqs. (24) and (3) indeed allows for a perfect fit of the data from the VES Dorofeev:2006fb () and BESIII Ablikim:2010kp () experiments. With data of even higher statistics, it might become possible to include the wave in the fit. In this case, one needs to determine the four subtraction constants of the DR representation in Eq. (3). A fifth subtraction constant would be introduced if we assumed a different highenergy behavior of : if a resonance with exotic quantum numbers coupling to exists (the search for which seems inconclusive so far Adolph:2014rpp (); Schott:2012wqa ()) and we assume the wave phase to approach instead of asymptotically, the wave would allow for a nonvanishing (constant) subtraction polynomial in Eq. (3), which cannot be removed by the transformation (3).
The inclusion of a wave contribution, which has been suggested in Ref. Escribano:2010wt () (in the form of the resonance) would require even higherorder subtraction polynomials.
4 Numerical solution of the dispersion relation
4.1 Iteration procedure
The numerical treatment of the integral equations (3) or (3) is a rather nontrivial matter, and specific care has to be taken in calculating the Omnès functions, the inhomogeneities with their complicated structure, as well as the dispersion integrals over singular integrands. All the details can be found in Ref. Schneider:2012nng ().
The solution of the integral equations is obtained by an iteration procedure: we start with arbitrary functions and , which we choose to be the respective Omnès functions; the final result is of course independent of the particular choice of these starting points. Then we calculate the inhomogeneities and , and insert them into the dispersion integrals (3) or (3). This procedure is repeated until sufficient convergence with respect to the input functions is reached. The iteration is observed to converge rather quickly after only a few steps.
The integral equations have a remarkable property that greatly reduces the numerical cost of the calculations: they are linear in the subtraction constants. Thus we can write (for the DR representation)
(28) 
where we have defined
(29) 
and analogously for the remaining basis functions. Each of the basis functions fulfills the decomposition Eq. (24), and we denote the corresponding singlevariable functions by , , etc., i.e.
(30) 
We can perform the iteration procedure separately for each of these while fixing the subtraction constants after the iteration.
4.2 Phase input
The crucial input in the dispersion relation consists of the and wave phase shifts and . Below the threshold of inelastic channels, the dispersion relation correctly describes rescattering effects according to Watson’s finalstate theorem with phases and that are equal to the phase shifts of elastic scattering. In principle, a coupledchannel analysis could be used to describe the process above the opening of inelastic channels, e.g. the explicit inclusion of intermediate states would provide a fully consistent description in the region of the and resonances; such a coupledchannel generalization of the Khuri–Treiman equations has recently been investigated for Albaladejo:2017hhj (). Alternatively, the singlechannel equations (3) remain valid if we promote and to effective phase shifts for this decay. As a full coupledchannel analysis is beyond the goal of this work, we construct effective phase shifts and quantify the uncertainties above the inelastic threshold. In such an effective onechannel problem, there are two extreme scenarios of the phase motion of at the resonance, depending on how strongly the system couples to strangeness Donoghue:1990xh (); Ananthanarayan:2004xy (): large strangeness production manifests itself as a peak at the position of the in the corresponding Omnès function, and thus the phase shift is increased by about while running through the resonance (this scenario is also realized in the elastic scattering phase shift). If the coupling to the channel with strangeness is weak, the corresponding Omnès function has a dip at the resonance position and the corresponding phase shift decreases (this is realized in the phase of the nonstrange scalar form factor of the pion). Scenarios in between these two extremes are conceivable.
For the input on the elastic phase shift, we use the results of very sophisticated analyses of the Roy (and similar) equations Caprini:2011ky (); GarciaMartin:2011cn (). As both analyses agree rather well, we only take one of these parametrizations Caprini:2011ky () into account. In Fig. 1, we show our phase , which agrees with the Roy solution Caprini:2011ky () below the inelastic threshold. The uncertainty due to the variation of the parameters in the Roy solution is shown as a light gray band labeled “lowenergy uncertainty.”
Now, the continuation into the inelastic region is modeled as follows. We calculate the waves for and in large PT at nexttoleading order (tree level) and unitarize this coupledchannel system with an Omnès matrix taken from Ref. Daub:2015xja (). The large PT representation depends on the lowenergy constants (LECs) and . We take their values from the most recent dispersive analysis of decays Colangelo:2015kha (),
(31) 
and vary each of them within its uncertainty. Adding the variations of the phase in quadrature generates the broad dark gray band labeled “highenergy uncertainty” in Fig. 1. This treatment correctly generates a smooth phase drop by with respect to the elastic scattering phase, and the uncertainty band covers a broad energy range for the position of this decrease. Still, the phase drops at sufficiently high energies such that the corresponding Omnès function, shown in Fig. 2, exhibits a peak at the position of the resonance. Asymptotically, we smoothly drive to a value of . We wish to emphasize that the role of the large input is not essential and only that of an auxiliary tool, which allows for a smooth construction that obeys the two desired features: the occurrence of an peak in accordance with expectations from scalarresonance models (see, e.g., Sect. 7.1), and an asymptotic phase of (as opposed to , say). The large highenergy uncertainty in Fig. 1 should safely cover a large variety of phases obeying these constraints. Note, finally, that in the physical region of the decay , the uncertainties of the phase and the Omnès function are small.
In the same spirit of an effective onechannel treatment, we consider isospinbreaking effects. In the isospin limit, our formalism applies identically to both the charged and the neutral processes and . In order to account for the most important isospinbreaking effects, we construct effective phase shifts for the neutral decay mode that have the correct thresholds and reproduce the expected nonanalytic cusp behavior, as we explain in detail in B.
For the phase shift , we take the phase of the scalar form factor of Ref. Albaladejo:2015aca () as input, shown in Fig. 3. In that reference, a – wave coupledchannel matrix has been constructed, to which chiral constraints Gasser:1984gg () have been imposed as well as experimental information on the and resonances. The remaining model uncertainty has been subsumed in the dependence on one single phase . The “central” solution corresponds to a parameter value of , while the uncertainty band is generated as an envelope of the solutions obtained by varying the parameter in the restricted range , which is compatible with chiral predictions for the scalar radius. The largest part of the highenergy uncertainty is generated by values of the parameter in the range , as shown in Fig. 3: while for all , the phase approaches its asymptotic value of very quickly above the resonance, this convergence becomes very slow and is extended over a vast energy range for . As this highenergy behavior turns out to affect the uncertainties in some (but not all) quantities extracted from data fits rather strongly, we will occasionally also refer to the reduced uncertainty, induced by the more restricted range . The Omnès function with an uncertainty band generated by the full variation is shown in Fig. 4. In particular, we observe a pronounced peak at the position of the resonance for all parameter values.
5 Determination of the subtraction constants
After having solved the integral equations numerically, we have to determine the free parameters in the dispersion relation, i.e. the subtraction constants , , and in the case of the DR representation, or , , , and in the case of DR. We summarize the experimental situation on Dalitz plots in Sect. 5.1. In Sect. 5.2, we discuss the results of fitting the subtraction constants to the most recent data sets.
5.1 Sampling of experimental Dalitz plots
In experimental analyses of the Dalitz plot, one defines symmetrized coordinates and according to
(32) 
where . The squared amplitude of the decay is then expanded in terms of these variables,
(33) 
and the parameters , , , are fitted to experimental data. Note that a nonzero value for the parameter (in ) would indicate violation of charge conjugation symmetry; there is no indication of a nonzero up to this point. In the following we consider two recent measurements of the charged final state . These determinations of the Dalitzplot parameters by the BESIII Ablikim:2010kp () and the VES Dorofeev:2006fb () collaborations currently feature the highest statistics. In Table 1 we have summarized some details and results of the two experiments.
BESIII Ablikim:2010kp ()  VES Dorofeev:2006fb ()  

# events  
# bins  
# bins 
For our analysis, we have generated pseudodata samples from the Dalitzplot distributions as measured by the two groups Kupscpriv (); the resulting Dalitzplot distributions are shown in Fig. 5. To check our results we have refitted the parametrization (33) to the synthesized data sets, and find agreement with the fit parameters of Table 1 within statistical uncertainties, as well as with the correlation matrix quoted in Ref. Ablikim:2010kp (). We note that the two data sets disagree on the parameter at the level; of course, it would be desirable that this experimental disagreement be resolved by future measurements.
5.2 Fitting experimental data
We proceed by fitting the subtraction constants, which are the free parameters in our dispersive representation of the amplitude, to the following data.

Dalitzplot distribution for the charged channel from BESIII Ablikim:2010kp () and VES Dorofeev:2006fb () experiments.

The partial decay width PDG ().
Note that we use real fit parameters: in principle the subtraction constants can have imaginary parts due to the complex discontinuity (3). However, since the imaginary parts of the subtraction constants are proportional to threeparticlecut contributions, and the available decay phase space of is small, the imaginary parts are so tiny that—given the precision of the data sets—their effect is entirely negligible (this is not the case for processes involving the decay of heavier mesons, compare e.g. Refs. Niecknig:2012sj (); Niecknig:2015ija ()).
In the following, we set up a scheme that allows us to fit both the experimental Dalitzplot distribution and the partial decay width simultaneously and avoids some strong correlations between the fitting parameters.^{4}^{4}4We write down formulae for the DR representation with subtraction constants , , and . The fitting procedure for the DR representation is completely analogous.
First, we perform the following transformation of the fit parameters:
(34) 
Hence, we write the squared amplitude as
(35) 
where
(36) 
and fix the arbitrary normalization of to
(37) 
This condition results in a quadratic equation for the rescaled subtraction constants , , and . We choose to express in terms of and . The experimental partial decay width now directly determines the normalization and has no influence on the parameters and