Dispersion-theoretical analysis of the D^+ -->K^- pi^+ pi^+ Dalitz plot

We study the Dalitz plot of the Cabibbo-favored charmed-meson decay $D^+\to K^-\pi^+\pi^+$ using dispersion theory. The formalism respects all constraints from analyticity and unitarity, and consistently describes final-state interactions between all three decay products. We employ pion-pion and pion-kaon phase shifts as input, and fit the pertinent subtraction constants to Dalitz plot data by the CLEO and FOCUS collaborations. Phase motions of resonant as well as nonresonant amplitudes are discussed, which should provide crucial input for future studies of CP violation in similar three-body charm decays.


Introduction
Heavy-flavor three-body decays into light mesons provide a valuable source for Standard Model tests and beyond. While they are driven, at short distances, by the weak interactions, their rich kinematic structure accessible in Dalitz plot distributions makes them a prime example for the application of modern tools of amplitude analysis [1]. A major motivation for the investigation of heavy-flavor decays is the study of CP violation, which manifests itself in the appearance of (weak) phases and requires the interference of different amplitudes with, at the same time, different phases in the strong final-state interactions (see, e.g., Ref. [2] for an in-depth overview). In contrast to (quasi-)two-body decays occurring at fixed total energies, three-body decays offer a resonance-rich environment with rapidly varying strong phases throughout the phase space available, which may strongly magnify the effects of CP violation in certain parts of the Dalitz plot.
Obviously, in order to turn the search for potentially very small CP -violating phases in such complicated hadronic environments into a precision instrument, it is inevitable to control the strong dynamics in the final state as accurately as possible, in a model-independent fashion that, however, incorporates a maximum of theoretical and phenomenological constraints. The traditional approach to model Dalitz plots in terms of isobars, i.e. a series of subsequent two-body decays, and describe the relevant line shapes in terms of Breit-Wigner (or Flatté) functions, has clear limitations: it fails to describe in particular the phase motion of the broad S-wave resonances such as the f 0 (500) in pion-pion or the K * 0 (800) in pion-kaon scattering (see e.g. Refs. [3,4] in the context of heavy-flavor decays), and neglects corrections beyond two-body rescattering in an unquantified manner.
It has therefore been advocated to employ the framework of dispersion theory for amplitude analyses [1], which is built on unitarity, maximal analyticity, and crossing symmetry. The dispersive framework adapted to study three-body decays was originally introduced by Khuri and Treiman for the decay K → 3π [5] and subsequently further developed [6][7][8][9][10]. The formalism has been resurrected in a modern form in Refs. [11,12]. The Khuri-Treiman equations are based on elastic two-body unitarity and explicitly generate crossed-channel rescattering between the three final-state particles. The equations are constructed by setting up dispersion relations for the crossed scattering processes, with a subsequent analytic continuation back into the decay region. This continuation is performed along the lines of the continuation of the perturbative triangle graph and is extensively discussed in Ref. [6].
Khuri-Treiman equations have been successfully applied to various low-energy meson decays, like e.g. η → 3π [11][12][13][14], ω/φ → 3π [15,16], or η ′ → ηππ [17]. In this work, we extend this formalism to three-body decays of open-charm mesons, analyzing the Cabibbofavored decays D + → K − π + π + /K 0 π 0 π + . As input we solely rely on ππ and πK phaseshift input. While these are not yet decay channels of major interest to study CP violation, the final-state interactions are going to be similar for others that are, such as the Cabibbosuppressed decays D → 3π/πKK. For the decays at hand, inelastic effects are small in large regions of the Dalitz plots, and therefore elastic unitarity provides a good approximation: the ππ channel allows for isospin 1 and 2 only, but no isoscalar components, which would necessitate a coupled-channel treatment, as a strong coupling to KK occurs. The major inelasticities in the πK channel are found to set in at the η ′ K threshold [18][19][20].
Thus with the high-statistics experimental data available [21][22][23], these decays provide a good test case to establish this dispersive framework in higher energy regions and set the path to Cabibbo-suppressed decays where traces of physics beyond the Standard Model may be searched for. Besides, it allows for a further test of low-energy πK and ππ dynamics as well as the importance of crossed-channel rescattering effects in three-body decays. It may also provide an insight into scattering phase shifts at higher energies in the future.
The decay under consideration has been the subject of a number of previous theoretical publications, focusing on different issues raised by the experimental results. One challenge is the proper treatment of the isospin 1/2 S-wave with the very broad, non-Breit-Wigner-shaped K * 0 (800) (or κ) resonance [24], and the inclusion of two scalar resonances K * 0 (800) and K * 0 (1430) in a way that conserves unitarity. Furthermore, the width of the K * 0 (1430) extracted from the experimental analyses in Refs. [21,22] is found to be inconsistent with PDG values [25]. In addition, the explicit comparison of the πK partial-wave phases extracted from these decays [23,26] with πK scattering results [27] seems to indicate deviations from Watson's final-state theorem.
Ref. [28] focuses on the isospin 1/2 S-wave final-state interactions, based on coupledchannel partial waves for Kπ, Kη, and Kη ′ constructed dispersively in Ref. [18]. Decay and scattering data could be reconciled, although no three-body rescattering effects, isospin 3/2 components, or ππ channel were included. Ref. [29] similarly observes mutual consistency of πK scattering and the D-meson decay, using related input to take two-body final-state interactions in the πK isospin 1/2 S-and P -wave into account in terms of the corresponding scalar and vector form factors. Furthermore, the short-distance weak interactions are described with the help of an effective Hamiltonian based on a factorization ansatz. Again, weak repulsive partial waves (of isospin 3/2 and in the ππ system) as well as crossedchannel rescattering are neglected. We mention that similar approaches, using dispersively constructed form factors for two-body rescattering, but neglecting third-particle interactions, have also been applied to B → Kππ decays [30,31].
In Ref. [32], a Faddeev-like equation is solved that builds up three-particle rescattering effects. The underlying two-particle πK amplitudes are obtained form unitarized chiral perturbation theory fitted to experimental data. The decay amplitude is simplified to include only the isospin 1/2 S-wave, aiming mainly at a study of the importance of rescattering effects and the reproduction of the experimental S-wave phases [23,26]. The model for the weak vertex has subsequently been improved [33]. Ref. [34] applies a similar approach with the addition of the isospin 3/2 πK S-wave, but is still restricted to S-waves only. The only theoretical analysis known to us with all relevant partial waves, three-particle rescattering effects, and effects of the intermediate stateK 0 π 0 π + included, is Ref. [35]. The author performs a full Dalitz plot analysis on pseudo data, which we will later compare to.
The outline of this article is as follows. Section 2 states some basic kinematical relations and shows both isospin and partial-wave decomposition of the decay amplitude in question. In Sec. 3, we derive the coupled dispersive integral equations and discuss how to solve these. Numerical results are shown in Sec. 4 and compared to experimental Dalitz plot data by the CLEO [21] and FOCUS [22] collaborations. We conclude our study in Sec. 5. Some technical details are relegated to the appendices.
We begin with the isospin and partial-wave decompositions of the decay amplitudes M −++ (D + → K − π + π + ) and M0 0+ (D + →K 0 π 0 π + ). We associate the isospin structure of the strong final-state current in Fig. 1 with the D + meson. Since oneūu/dd pair is strongly produced, the associated isospin of the D meson is given by I = 3/2, I z = 3/2. Thus the isospin decomposition of the respective (crossed) scattering processes reads where F I denotes the amplitude with definite isospin I. The decay amplitudes are given by We can write down a symmetrized partial-wave expansion simultaneously in s-, t-, and u-channels (the precise relation of which to proper partial waves in a single channel will be discussed below). With the expansion in partial-wave amplitudes truncated at the D-wave The gray vertex stands for the crossed decay amplitude D + π − → K i π j denoted by M ij+ and the white vertex the K i π j → K − π + scattering amplitude denoted by T ij,−+ . The dashed line gives the contribution to the discontinuity [41]. The other channels follow analogously.
for πK final states and the P -wave for ππ, we obtain where the single-variable amplitudes F I L have definite isospin I and angular momentum L in the channel associated with the Mandelstam variable featuring as their argument. Note that the inclusion of D-waves is somewhat heuristic: in order to rigorously prove the symmetrized decomposition (2.5) in the spirit of the so-called reconstruction theorem [12,[36][37][38][39][40], one needs to include a subtraction polynomial of higher order (i.e., a larger number of unknown parameters) than what we will allow for below. We mainly want to retain the πK D-wave to test the effect of the K * 2 (1430) resonance, which is kinematically accessible in the decay phase space. The way we implement this approximately will be discussed in Sec. 3.3.

Unitarity and Omnès solution
We begin with the dispersive treatment of the associated scattering processes linked to the decay by crossing symmetry, D +π → Kπ and D +K → ππ. The D-meson mass is artificially set to M D < M K + 2M π such that the corresponding decay is kinematically forbidden. The simpler analytic structure of these scattering processes can be exploited to construct dispersion relations for the single-variable amplitudes valid for s, t > (M D + M π ) 2 and u > (M D + M K ) 2 , respectively. The analytic continuation back to the physical Dmeson mass as well as into the kinematic region yields the anticipated decay amplitudes [6]. We demonstrate the framework for the example of the s-channel processes; the t-and uchannel amplitudes are constructed analogously. Elastic unitarity gives for the discontinuity (see Fig. 2  where the sum runs over isospin and angular momentum components I and L. Furthermore we use the Clebsch-Gordan coefficients a I,L , Legendre polynomials P L (z), and the corresponding partial waves t I L (s) and f I L (s). 1 Exploiting the unitarity relation for elastic πK and ππ scattering we obtain the following partial-wave unitarity relations 1 Note that in contrast to the definition of the single-variable amplitudes in Eq. (2.5), we have not defined the partial waves in Eq. (3.2) to be free of kinematical zeros. This is independent of the singularities these partial waves display at the corresponding pseudo-thresholds or upper limits of the physical decay region, s = (MD − Mπ) 2 or u = (MD − MK ) 2 , which are well understood, see e.g. Ref. [15] or the discussion in Ref. [42] in a perturbative context.
where δ I L (s) denotes the elastic final-state scattering phase shift. The thresholds in the different channels are s th = t th = (M K + M π ) 2 for πK and u th = 4M 2 π for ππ scattering, respectively. Since the discontinuity of f I L and the according single-variable amplitude κ L F I L coincide on the right-hand cut, we have where we have introduced the inhomogeneitiesF I L (s) that are free of discontinuities on the right-hand cut by construction. They incorporate the left-hand cut contributions and will be further discussed in Sec The inhomogeneous solution is obtained by a product ansatz where P I L (s) is now a polynomial of order n − 1, and the number of subtractions n is chosen such that the convergence of the dispersion integral is guaranteed.
As our approach relies on elastic unitarity (see Ref. [44] for a generalization of the Khuri-Treiman formalism to coupled channels), the formalism breaks down when inelastic channels become important. We assume that Watson's theorem [45] is a good approximation up to the η ′ K threshold in the πK channel. Inelastic effects in the prominent πK S-wave systems are found to become sizable above the η ′ K threshold [18][19][20]. The main inelastic contributions in the isospin 1/2 P -wave come from the πK * and ρK channels, which become noticeable in the energy region where they couple to K * (1410) and K * (1690) [20]. In all exotic partial waves, i.e. the isospin 2 ππ system as well as I = 3/2 πK partial waves, inelastic effects are assumed to be negligible.

Inhomogeneities
With the scattering phase shifts given as fixed input, the only quantities left in the dispersion integrals Eqs. (3.7) are the inhomogeneitiesF I L , which are determined as the projections of the crossed-channel amplitudes onto the considered channel. They can be re-expressed in terms of the single-variable amplitudes F I L (x), such that we obtain integral equations that can be solved for the F I L (x). With the aid of Eq. (3.4) we find where M Ix ijk (x, z x ) denotes the projection of the full decay amplitude M ijk (x, z x ) onto isospin I x eigenfunctions in the x-channel. One term of the projection integral over , such that the right-hand-cut discontinuity is canceled. The inhomogeneities are then indeed free of discontinuities on the right-hand cut as anticipated. The resulting inhomogeneities are given in Appendix A.
The interpretation of Eq. (3.8) as an angular integration is valid in the scattering region and needs to be analytically continued into the unphysical and decay regions. Performing the angular integration naively in the decay region results in crossing the unitarity cut. The prescription on how to perform the continuation has been extensively discussed in Ref. [6], motivated by the continuation of the (perturbative) triangle graph into the decay region. It ultimately leads to the prescription M 2 D → M 2 D + iǫ, which allows one to derive an integration path that avoids the unitarity cut.

Number of subtraction constants
The minimal number of subtractions needed is dictated by the asymptotic behavior of the integrands in Eqs. (3.7). The decay amplitude and thus the inhomogeneities are assumed to grow at most linearly asymptotically, loosely based on the Froissart bound [46]. Assuming the phase shifts to approach constant values δ I L (∞) for large energies, the Omnès functions With the following assumption for the phase shifts δ I L at high energies: we need two subtractions for F 2 0 , F 1 1 , and F 3/2 0 , four subtractions for F 1/2 0 , and one subtraction for F 1/2 1 to obtain convergent dispersion integrals. Note that the difference in the number of subtractions for F 1 1 and F 1/2 1 , despite identical phase asymptotics, is due to the different kinematic prefactors for P -waves with equal and unequal masses, see Eq. (2.5). F 3/2 1 needs no subtraction, but as the πK isospin 3/2 P -wave phase shift is very small and assumed to vanish at high energies, we neglect it altogether. Similarly, also the I = 3/2 D-wave is put to zero.
The inclusion of the D-wave F 1/2 2 is delicate. Formally it requires no subtractions, but the kinematical pre-function corresponding to the L = 2 Legendre polynomial, multiplied with the required momentum factors to make it free of kinematical singularities, see Eq. (2.5), violates the assumed high-energy behavior of the decay amplitude and thus of all inhomogeneities. Therefore we will follow a "hybrid approach" for the D-wave: we will only consider the projections of S-and P -waves of other channels in order to generate the D-wave inhomogeneity, but will exclude D-wave projections, thus eschewing the need for further subtractions. This is loosely motivated by analogous observations in low-energy processes calculated in chiral perturbation theory, where higher partial waves are dominated by crossed-channel loop diagrams that correspond to low partial waves in those crossed channels.
In total we have eleven subtraction constants. However, the resulting representations of the decay amplitudes Eqs. (2.5) are not unique due to the linear dependence of the Mandelstam variables s, t, and u: One can construct polynomial contributions to the single-variable amplitudes that leave the complete decay amplitudes M −++ (s, t, u) and M0 0+ (s, t, u) invariant; this is obvious in a standard dispersive representation, however slightly less trivial to demonstrate in the Omnès representations discussed above [47]. The polynomial coefficients can be tuned such that a maximal number of subtraction constants is eliminated to obtain a linearly independent set. These polynomials span the so-called invariance group of the decay amplitudes. Details are discussed in Appendix B. We choose to eliminate the subtraction constants in the nonresonant I = 3/2 πK and I = 2 ππ Swaves, the rationale being solely to retain them in presumably large, resonant partial waves. This leaves seven linearly independent complex subtraction constants, . (3.10) The subtraction constants cannot be determined in the framework of dispersion theory and have to be obtained either by matching to a more fundamental dynamical theory, or, as in this work, by a fit to experimental data. The solution space of the coupled system Eq. (3.10) has thus dimension seven, corresponding to the seven complex subtraction constants. Since the equations depend linearly on the subtraction constants, it is convenient to choose seven independent basis sets and solve the equations for each of these sets. We call those solutions basis functions. In particular, we choose for the ith basis function M i (s, t, u) the set of subtraction constants The basis functions are entirely determined by the phase shift input, as well as the masses of all particles involved (taken from Ref. [25]). The phase shifts are obtained from solutions of the ππ Roy equations by both the Bern [48][49][50] and the Madrid-Kraków [51] groups, as well as the Roy-Steiner equations for πK scattering solved by the Orsay group [52].

Solution strategy
In this section we discuss different solution strategies of the Khuri-Treiman-type equations (3.10), their issues, and present our new solution strategy.
The standard solution strategy for the linear coupled double integral equations (3.10) has been an iteration procedure as performed for example in Refs. [13,15,17] or numerically faster with the introduction of integral kernels in Ref. [53]. Starting with an arbitrary input for the single-variable amplitudes (e.g. just the Omnès functions), the inhomogeneities are evaluated; with these the dispersion integrals are determined to obtain a new set of single-variable amplitudes. This cycle is repeated until satisfactory convergence is reached. Unfortunately, the convergence of this iterative procedure is not always guaranteed, depending on the mass of the decaying particle and the number of subtractions: for larger decay masses and more subtractions applied, the corrections in each iteration step can be too large to reach the fixed-point solution. We find this to be the case in the D-meson decays considered here.
This necessitates a different solution strategy. Since the set of integral equations is linear in the single-variable amplitudes it is convenient to set it up in the form of a matrix equation instead. Provided that the matrix is invertible a unique solution exists. One such inversion strategy is known as the Pasquier inversion [9] (see Ref. [54] for a recent comparison of Pasquier inversion and iterative solution), where a method to reduce the double integral equation to a single integral equation is introduced. The procedure involves the deformation of the integral contours of both integrals, allowing one to interchange the order of integrations such that a unique kernel function is obtained. The coupled single integral equations thus obtained do allow for a direct solution via matrix inversion.
We will follow a slightly modified strategy, constructing a matrix equation without performing a Pasquier inversion. In this context it is beneficial to solve for the inhomogeneities instead of the single-variable amplitudes, the advantage being that the inhomogeneities need to be evaluated only on the right-hand-cuts. The single-variable amplitudes themselves can be obtained in the whole complex plain in a straightforward manner by performing the dispersion integral over the inhomogeneities once.
To illustrate the solution strategy we limit ourselves to one hypothetical inhomogeneity equation without any loss of generality, and focus on the functionsF(s) ≡F L (s)κ 2L+1 (s) that are free of singularities at the pseudo-threshold or upper limit of the kinematically accessible decay region (which is a zero in κ(s)). Inserting Eq. (3.7) into Eq. (3.12) yields The function A(s) contains the dependence on the subtraction polynomial, while the integration kernel K(s, x) is independent of any subtraction constants: (3.14) Equation (3.13) is thus a linear integral equation forF (s), to be solved for a given set of subtraction constants. Discretizing Eq. (3.13) yields which is solved by matrix inversion; the numerical treatment is relegated to Appendix C.

Numerical results and experimental comparison
Solving the coupled integral Eqs. The error bands in Figs. 3 and 4 are determined by a conservative error estimate of the phase shifts: For the S-wave πK and ππ phases the error is assumed to rise linearly from zero at threshold to ±20 • at 2 GeV. Beyond 2 GeV the error is fixed to ±20 • . The πK isospin 1/2 P -and D-wave phase errors and ππ P -wave phase errors are similarly obtained, with the only difference that the linear rise of the error sets in after the K * (892), K * 2 (1430), and ρ(770) resonances, respectively. In the ππ P -wave case we additionally vary between the phase-shift data from Refs. [48][49][50][51].
In the following we compare our theoretical decay amplitude to the experimental D + → K − π + π + Dalitz plot data from the CLEO [21] and FOCUS [22] collaborations. Exploiting the symmetry of the process under the interchange of the two pions, we can restrict the comparison to the region s < t by mirroring the remaining half of the Dalitz plot into this region.   distribution function was used for the fit analogously to the experimental analyses with (s i , t i ) being the center of the corresponding bin and 2δ the bin width, ǫ(s, t) the efficiency parametrization, B(s, t) the background parametrization, N sig and N B normalization constants such that the background and signal term are normalized to unity, and the signal fraction f sig .
We minimize the following χ 2 , where N is the number of events, the sum runs over the number of bins and the error on the binned data is assumed to be purely statistical. In addition to the full dispersive representation Eq. (3.10), we also fit a simplified decay amplitude to data, which is given by a sum of Omnès functions multiplied by polynomials: where the c ′ i are again complex fit constants. Equation (4.3) emulates a dispersively improved isobar model that neglects any crossed-channel rescattering effects. The number of polynomial fit constants is chosen to resemble the number of degrees of freedom in the full dispersive result Eq. (3.10) as far as possible; with certain caveats that preclude an immediate quantification of three-particle rescattering effects in the same straightforward way as performed for φ → 3π decays in Ref. [15]. In Eq. (3.10), two subtraction constants c 0 and c 1 are contained in the ππ P -wave, which only contributes indirectly via the intermediate statē K 0 π 0 π + to the decay and thus does not show up in the pure Omnès amplitude Eq. (4.3). In addition, every Omnès function in Eq. (4.3) needs at least a normalization constant to adjust the strength of individual amplitudes, while some single-variable amplitudes do not have any subtraction constants. Finally, once the D-wave is included we have one additional complex fit parameter c ′ 7 in the pure Omnès fits. For that reason we consider both Omnès and the full dispersive fits without (Omnès 1, full 1) and with D-wave (Omnès 2, full 2).
We have the freedom to fix one subtraction constant, as both the overall normalization and the overall phase are arbitrary and factorized out; we choose c 2 = c ′ 2 = 1. This leaves 13 (15) real fit constants for the full / Omnès fits.
Following experimental custom, we will employ so-called fit fractions to characterize the relative importance of various single-variable functions. These are defined in the following way where the P J denote the angular prefactors of the corresponding single-variable amplitudes in the total amplitude. The integration runs over the fitted Dalitz plot region. In general these fit fractions are not unique due to the freedom of adding an element of the invariance group Eq. (B.1); the projections onto partial-wave amplitudes then will lead to different fit fractions.

Comparison to the CLEO data
The Dalitz plot measured by the CLEO collaboration [21] contains 140793 events. The efficiency and background parametrizations are given explicitly. 2 Our fit results are summarized in Table 1, together with the fit fractions in Table 2. In the full dispersive fits (full fits 1/2), the resulting values for the subtraction constants in Table 1 have similar order of magnitude with the exception c 6 , which is rather small. This can be understood by the large F 1/2 1 single-variable amplitude in this particular basis function (see Fig. 4). Further- more the phases of the F 1 1 subtraction constants (c 0 , c 1 ) nearly agree modulo π. The same holds for the F 1/2 0 subtraction constants (c 2 to c 5 ) especially for the full fit 2. This suggests that with overall phases factorized, the subtraction constants for the F 1 1 and likewise the F 1/2 0 amplitude are almost real. The differences of the single-variable amplitude phases to the elastic phase shifts depicted in Fig. 7 are thus predominantly due to the dispersion integrals, i.e. the crossed-channel rescattering effects.
Including the D-wave improves the χ 2 /d.o.f. slightly from 1.18 ± 0.03 to 1.10 ± 0.02. Note that in the full dispersive representation, no additional fit constants are introduced when the D-wave is added. The inclusion of the D-wave does not change the phases of most subtraction constants beyond their uncertainties, with the exception of c 6 ; the magnitudes, in contrast, change significantly for almost all subtractions. Considering the fit fractions in Table 2, we observe that the inclusion of the D-wave in the full fit 2 reduces the highly destructive interference between the two S-wave amplitudes in the πK channel. We wish to point out that also in Ref. [22], a large cancellation between the isospin 1/2 and isospin 3/2 S-wave components of −164% is seen, with individual fit fractions of (207 ± 24)% and (40 ± 9)%, respectively, which show a comparable behavior to our full fit 1. Although the fit fraction of the D-wave itself is very small, it thus has a rather large impact on the Sand P -waves. A similar phenomenon is seen in Ref. [21] where the fit quality deteriorates considerably when removing the small D-wave. Although we do not fit the whole Dalitz plot, the fit fractions for the resonant single-variable amplitudes for F 1/2 0 , F 1/2 1 and F 1/2 2 agree well with the results from Refs. [21,22]. The F 2 0 fit fraction corresponds to the isospin 2 ππ S-wave component of FF ≈ (9.8 . . . 15.5)% found in Ref. [21] within different fit models, and together with the fit fraction of F 3/2 0 agrees with the nonresonant contribution found in Ref. [22] of FF ≈ (29.7 ± 4.5)%.
Although the Omnès fits (Omnès 1, 2) yield overall similar χ 2 results, the strengths of the individual amplitudes shown in Table 2 are highly implausible and probably sufficient to reject this model. In particular the contribution of the nonresonant isospin 3/2 πK S-wave is vastly beyond all reasonable expectations, and cannot be justified. In contrast to the full fit, this situation is not ameliorated significantly by including the D-wave. We conclude that crossed-channel rescattering effects are essential to obtain sensible fit fractions.
1.20 ± 0.01 1.21 ± 0.02 1.25 ± 0.02 1.17 ± 0.01 Table 3. Fit to FOCUS data: Numerical fit results for the subtraction constants c i and c ′ i and the corresponding χ 2 /d.o.f.. The same four fit scenarios as in Table 1 are considered. The errors on the parameters are evaluated by varying the basis functions within their error bands.
The resulting Dalitz plot as well as a one-dimensional representation in terms of slices through it are displayed in Fig. 5. The bin numbering for the latter is organized in terms of t-slices for constant s, subsequently glued together with the next slice of higher s. We evaluate the event distribution function Eq. (4.1) over each bin and compare to experimental data. The rather small error band on the fit results suggests that the uncertainty in the basis functions is largely compensated by interference effects between the different single-variable amplitudes, as well as by corresponding variations in the fitted subtraction constants.

Comparison to the FOCUS data
The FOCUS Dalitz plot data [22] includes 52460 ± 245 signal and 1897 ± 39 background events. With the resulting signal fraction of ∼ 96.5% we perform the full and Omnès fits as above. Table 3 summarizes the fit results together with the fit fractions in Table 4. The overall picture is very similar to the CLEO fit results with a slightly bigger χ 2 /d.o.f. ≈ 1.2. The Omnès fits again result in nonphysical fit fractions (see Table 4), and from here on we will only compare the full fits of both experimental data sets. Starting with the fit without D-wave (full fit 1) we observe similar moduli of the subtraction constants compared to the CLEO results, however the phases do differ. The fit does not show the large destructive interference effects between the isospin 1/2 and isospin 3/2 S-wave that we find in the CLEO fit.
No improvement in the χ 2 /d.o.f. is observed when we include the D-wave (full fit 2). However the contribution from the nonresonant amplitudes, the isospin 2 and isospin 3/2 S-waves, are reduced (see Table 4). The fit fractions of the full fit 2 differ slightly from the CLEO fits; in particular the nonresonant S-waves contribute less in the FOCUS data.
In the full fit 2 the phases of the F 1 1 subtraction constants persist to nearly agree modulo π; the same holds for F 1/2 0 subtraction constants. It is reassuring that the overall picture of the phases of various subtraction constants is consistent in the full fit 2 results for both CLEO and FOCUS.
In Fig. 7, we compare moduli and phases of the resulting single-variable amplitudes as fitted to the two data sets; the phases are also compared to the input phase shifts used in the Omnès functions. The resulting phase motions largely agree in the two analyses within uncertainties, with the possible exception of some deviations in F 1/2 0 in the region of the K * 0 (800) resonance, where the phase extracted from the CLEO fit rises more quickly. There are significant deviations from the input phase shifts throughout: there is no naive realization of Watson's theorem in the presence of three-body rescattering effects, see e.g. recent discussions in Refs. [54,56]. This is also the explanation for the observed discrepancy of the πK I = 1/2 S-wave phase as extracted from these decays by the E791 [23] and FOCUS [26] collaborations, compared to the scattering phase-shift analyses [27]: while the phase shift rises to about 67 • -97 • at √ s = 1.3 GeV [52], the experimental analyses of Ddecay data suggest an increase in the phase from threshold by about 133 • -164 • (read off via Ref. [32]). Figure 7 shows that in the dispersive formalism, the phase at 1.3 GeV is about 182 • -198 • (CLEO) or 170 • -183 • (FOCUS)-even larger than found in Refs. [23,26]. We emphasize that these results are based on a formalism that uses the scattering phase shifts [52] as input: the deviations in the decay amplitude S-wave are due to complex phases  . The overall normalization is chosen such that the absolute values in the K * (892) peak agree. Right column: Phases of the single-variable amplitudes (CLEO: red, FOCUS: blue) and input scattering phases (black) in radiant. The phases are fixed to zero at the two-particle (ππ, πK) thresholds. The dotted lines visualize the fitted area; for the πK amplitudes from threshold to the η ′ K threshold and the full phase space for the ππ amplitudes.
induced by three-body rescattering effects.
In general, the corrections compared to input phase shifts are smallest for narrow resonances, in particular in the I = 1/2 πK P -and D-waves. The largest phase differences are observed in the nonresonant amplitudes, where the phases of F 2 0 and F

3/2 0
show a 2π rise due to zeros in imaginary or real parts close to threshold in individual basis functions. Note how these seemingly drastic differences are accompanied by very small absolute magnitudes of the amplitudes in question: in view of the aim to control the phase behavior of the complete, combined decay amplitude accurately, these specific deviations are still rather small.
Turning to the moduli of the single-variable amplitudes, the relative strength of the K * (892) resonance (in F 1/2 1 ) compared with the K * 2 (1430) (in F 1/2 2 ) agrees between CLEO and FOCUS fits. However the dip in the F 1/2 0 amplitude is shifted to higher energies in the FOCUS fit and slightly more pronounced. The moduli of the nonresonant amplitudes F 3/2 0 and F 2 0 turn out to be smaller in the FOCUS fit, which is also underlined by the fit fractions (compare Tables 2 and 4).

Comparison to other approaches
So far the only theoretical approach known to us that includes all relevant partial waves, three-particle rescattering effects, and the isospin coupled intermediate stateK 0 π 0 π + is Ref. [35]. The treatment is based on a unitary coupled-channel framework. The two-particle rescattering contributions are fixed by the πK and ππ scattering data, phases, and moduli. Three-body rescattering effects are generated by solving a Faddeev equation. In addition to the three-body rescattering a three-body potential, based on hidden local symmetry, is introduced modeling vector meson exchanges. The author studies the influence of individual rescattering contributions by considering different fit scenarios; crossed-channel rescattering effects and three-body potential turned off (isobar fit), three-body potential turned off (Z fit), and the full fit. An additional contact term breaking unitarity is allowed for, which in the full fit turns out to be negligible. The decay amplitude depends on 27 to 39 degrees of freedom depending on the considered fit model, which is more than twice the number of parameters included in our full fit.
To compare the fit fractions obtained in Ref. [35], we note that the isobar fit theoretically compares closest to our Omnès fits, while the Z fit does to our full fits. However the isobar fit has a large contribution from the unitarity-breaking contact term (considered as a "background" contribution) of 17.7%, such that a direct comparison is not sensible. Concerning the full and the Z fit, a large destructive interference between the isospin 1/2 and isospin 3/2 S-waves is seen, similar to our CLEO fit 1 configuration. The isospin 1/2 P -waves are of similar size, ∼ 15% compared to our 10 − 14%, but the ππ S-wave contribution is smaller (1.8 − 3.8%) than our contributions in either full fit 1 or CLEO full fit 2. It agrees only with the FOCUS full fit 2. -Concerning this comparison, we should stress once more that in contrast to Ref. [35], we do not fit the full Dalitz plot.
Unfortunately the improvement due to crossed-channel rescattering cannot be quantified in a simple way in Ref. [35] either. The improvement going from the isobar to the Z and then further to the full model can also be due to the introduction of further degrees of freedom; as discussed above, we encounter a similar problem in our analysis. However the background term, which gives an indication for missing physics, reduces dramatically once the crossed-channel rescattering effects and the coupled intermediate stateK 0 π 0 π + are included. This is a similar conclusion as drawn from the dispersive analysis of φ → 3π Dalitz plots [15], which rendered phenomenological contact terms [57,58] superfluous.

Conclusion
In this paper we have analyzed the D + → K − π + π + decay with a dispersive framework based on the Khuri-Treiman formalism that satisfies analyticity, unitarity, crossing symmetry, and includes crossed-channel rescattering among the three final-state particles.
The theoretical decay amplitude depends on seven complex subtraction constants, one of which can be absorbed into overall phase and normalization of the amplitude. The remaining parameters are fitted to the experimental Dalitz plot data from the CLEO [21] and FOCUS [22] collaborations, restricting the kinematic region to below the η ′ K threshold, where the elastic approximation is assumed to work well. We have considered different fit scenarios with (full) and without crossed-channel rescattering effects (Omnès), as well as with and without the πK isospin 1/2 D-wave. Although the Omnès fits give reasonable χ 2 /d.o.f., we obtain large destructive interferences between single-variable amplitudes, which manifest themselves in unphysical fit fractions. The full fits result in good χ 2 /d.o.f. around 1.1 for the CLEO data (1.2 for the FOCUS data), with sensible fit fractions throughout. Including the πK isospin 1/2 D-wave does not significantly improve the χ 2 /d.o.f., however the fit fractions of the nonresonant waves are reduced, giving small interference effects between the single-variable amplitudes. We have shown that we can describe the D + → K − π + π + Dalitz plot data in the region where we deem elastic unitarity to hold approximately, solely relying on ππ and πK scattering phase shift input and exploiting the constraints of dispersion theory.
Three-body rescattering effects suspends any strict relation between the phase of the decay partial waves and scattering phase shifts: we have shown that the significantly stronger rise of the πK S-wave phases, as observed in analyses of these D-meson decays [23,26] in comparison to phase shift data, can be understood at least qualitatively in the framework of Khuri-Treiman equations.
We have simultaneously constructed the formalism for the decay D + →K 0 π 0 π + , which is directly related to D + → K − π + π + by charge exchange and can be constructed from different linear combinations of the same (isospin) amplitudes. This second decay channel has recently been measured by the BESIII collaboration [59]. A simultaneous analysis of both Dalitz plots will further exploit the predictive power of the dispersive formalism; due to the direct contribution of the ππ P -wave in the π 0 π + (as opposed to the π + π + ) final state, we expect to find stronger constraints on the subtraction constants featuring directly in the corresponding amplitude. The pertinent investigation is in progress [60].

A Inhomogeneities
In this appendix the inhomogeneities are calculated from Eq. (3.8). To demonstrate the procedure we will perform the calculation explicitly in the case of f I=1/2 L (s), We start with the projection of the decay amplitude M −++ , Eq. (2.5), onto isospin eigenstates in the s-channel. We introduce the following crossing matrices: and so on, where M I x is the isospin I eigenstate in the x-channel and X II ′ xy the crossing matrix for the transition from channel y to x, where I and I ′ are the matrix component indices. We obtain the following explicit forms: The t-channel and u-channel single-variable amplitudes can be split, with the aid of the crossing matrices, into I s = 1/2 and I s = 3/2 contributions,  We immediately obtain z n f us = z n f ut , z n f ts = z n f st , and z n f tu = (−1) n z n f su .
The angular average integration is straightforwardly performed in the scattering region. The continuation to the decay region, where the naive integration would cross the righthand cut, has been discussed extensively before [6,15]. We now perform the partial-wave projection with the Legendre polynomials P L (z s ). For the S-wave we obtain x ∈ {s, t}. Thus from Eq. (A.9), the inhomogeneity can be immediately read off from the relation √ 3M   where in addition to Eq. (A.10) we have used (A.12)

B Invariance group matching
In this appendix, we study the polynomial ambiguities in the decomposition of the total decay amplitudes Eq. (2.5) into single-variable functions, dubbed "invariance group". We wish to determine the polynomial at most linear in the Mandelstam variables that can be added to the different single-variable amplitudes, leaving the total decay amplitudes Eq. (2.5) invariant. For this purpose, we make use of the relation s + t + u = 3s 0 = M 2 D + M 2 K + 2M 2 π . It is easy to check that adding the following terms to the various S-waves as well as the ππ P -wave: As an example we will study the single-variable amplitude F 2 0 with F 2 inv 0 (s) = a 0 + b 0 s. We obtain π 2 0 (s) = a 0 + b 0 + a 0 (ω 2 0 ) 1 s + (ω 2 0 ) 1 − 1 π Using the sum rule value for (ω 2 0 ) 1 we find π 2 0 (s) = a 0 + b 0 + a 0 (ω 2 0 ) 1 s . (B.10) The other subtraction polynomials are obtained in an analogous way and read π 2 0 (s) = a 0 + (b 0 + a 0 (ω 2 0 ) 1 )s , π 1 1 (s) = − with no contributions to the πK P -and D-waves. Polynomial terms with higher order than the subtraction polynomials of the corresponding amplitudes (see Sec. 3.3) have been omitted.