Consistent Dalitz plot analysis of Cabibbo-favored $D^+ \to \bar{K} \pi \pi^+$ decays

We resume the study of the Cabibbo-favored charmed-meson decays $D^+ \to \bar{K} \pi \pi^+$ in a dispersive framework that satisfies unitarity, analyticity, and crossing symmetry by construction. The formalism explicitly describes the strong final-state interactions between all three decay products and relies on pion-pion and pion-kaon phase shift input. For the first time, we show that the $D^+ \to K_S \pi^0 \pi^+$ Dalitz plot obtained by the BESIII collaboration as well as the $D^+ \to K^- \pi^+ \pi^+$ Dalitz plot data by CLEO and FOCUS can be described consistently, exploiting the isospin relation between the two coupled decay channels that provides better constraints on the subtraction constants.


Introduction
Three-body decays of heavy mesons provide a powerful mean for Standard Model tests and beyond. Due to their richer kinematic structure as compared to two-body decays, they have a wide range of applications e.g. in hadron spectroscopy, studies of hadronic final-state interactions, or CP-violation studies. For these investigations, a thorough understanding of the strong final-state interactions is mandatory, necessitating the methods of amplitude analysis [1].
One particular issue in this area are sensible parametrizations of scalar partial waves. It is well known that the lowest-lying scalars cannot be described in terms of Breit-Wigner functions. In particular the lightest pion-pion and pion-kaon scalar resonances, the f 0 (500) (or σ) and the K * 0 (800) (or κ), are associated with poles in the complex energy plane too far away from the real axis to allow for any such simplistic description [2][3][4][5]. This has been recognized also in the context of heavy-meson decays, where more appropriate descriptions in terms of scalar form factors have been applied [6][7][8], using model-independent methods such as chiral perturbation theory and dispersive techniques.
A different aspect of three-body decays that requires refinement compared to a Breit-Wigner type isobar description lies precisely in the presence of a third hadron in the final state: crossed-channel rescattering effects will necessarily modify the spectral forms of different resonances, and the extent to which this is the case (and can potentially be regarded as a "small correction") is too seldomly investigated explicitly. One wellestablished theoretical tool to study such modifications are the Khuri-Treiman equations [9][10][11][12][13][14] (see also the recent lecture notes in Ref. [15]), originally invented to study K → 3π decays in a manner consistent with the constraints from analyticity (i.e., causality) and unitarity (i.e., probability conservation). They were resurrected and applied extensively to study η → 3π [16][17][18][19][20][21], and have subsequently also been applied to other low-energy three-body decays such as ω/φ → 3π [22,23] or η ′ → ηππ [24,25].
Recently, we have applied Khuri-Treiman equations to the Cabibbo-favored charmed-meson decays D + → K − π + π + /K 0 π 0 π + [26]. In comparison to the light-meson decays mentioned above, the corresponding Dalitz plot is significantly larger, with a much richer spectrum of partial waves/resonances contributing. In particular the decay D + → K − π + π + has been theoretically studied frequently before [27][28][29][30][31] (see also Refs. [32,33] for analyses of similar B-meson decays), using various approximations in the description of final-state interactions. This corresponds to the rather good data situation for that channel [34][35][36][37]. However, a consistent, combined investigation of both final states is highly desirable, as they are coupled to each other by simple charge-exchange rescattering, but only the partially neutral final state allows for the observation of resonances in the pion-pion system (essentially the ρ(770)), while the π + π + system is necessarily a nonresonant isospin I = 2 state. In simple isobar models as conventionally used in experimental analyses, neglecting the interaction of all three final-state particles, the relation between both channels is therefore obviously lost. In addition to Ref. [26], such a combined theoretical analysis has only been performed in Ref. [38], in the latter case using Faddeev equations to generate threebody rescattering effects. With the advent of experimental data on D + →K 0 π 0 π + , courtesy of the BESIII collaboration [39], we are now in the position for the first time to test our theoretical approach for consistency, using real data for both channels. This is the purpose of the current letter. We briefly summarize the dispersion-theoretical formalism developed in Ref. [26] in Sect. 2, before performing fits to the new BESIII data as well as combined fits for both final states in Sect. 3. We find the need to somewhat improve on the amplitude representation in partic-ular with regards to the D-wave, which is discussed in Sect. 4. We conclude our study in Sect. 5.

Kinematics, decay amplitude, dispersive representation
We define the Mandelstam variables of the three-particle decays The corresponding crossed-channel scattering angles are given by where the Källén function is defined by λ(x, y, z) ≡ x 2 + y 2 + z 2 − 2(xy + xz + yz). The decay amplitudes are decomposed into amplitudes depending on one Mandelstam variable only along the lines of the so-called reconstruction theorem [17,[40][41][42][43][44]. The explicit decompositions for the decay channels considered here have been performed in Ref. [26] and read for the D + →K 0 π 0 π + decay and for D + → K − π + π + . The strong final-state interactions of both decay channels are isospin-related and can therefore be described by the same single-variable amplitudes F I L of isospin I and angular momentum L. The above decomposition is consistently truncated beyond D-waves, and we have neglected exotic, nonresonant partial waves beyond the S -waves, i.e. the tiny πK P-and D-wave amplitudes of isospin I = 3/2 as well as the ππ I = 2 D-wave. 1 The F I L satisfy the following elastic unitarity relations: where x th denotes the elastic threshold of the channel considered, and δ I L (x) the corresponding ππ/πK scattering phase shift input of isospin I and angular momentum L taken from Refs. [45][46][47][48] (see also Ref. [49] for a new analysis of πK scattering). We attribute uncertainty bands to all phase shifts that rise linearly from zero at threshold (for the S -waves) or the position of the first resonance (for the P-and D-waves) to ±20 • at 2 GeV; see Fig. 2 (right column) below. The inhomogeneitieŝ F I L (x) are given by the subsequent partial-wave projectionŝ and give rise to crossed-channel rescattering contributions. They are calculated explicitly for all channels in Ref. [26]. The constants a i jk I,L denote the Clebsch-Gordan coefficients corresponding to the single-variable amplitudes for the final states i jk ∈ {00+, − + +}. The solutions to these unitarity relations, Eq. (6), are given in the form of inhomogeneous Omnès solutions: , with Ω I L (x) the corresponding Omnès functions While built on the requirement of fulfilling two-body unitarity, Eq. (6), the amplitude representation remarkably also fulfills the constraints of three-body unitarity [12,15] (compare also Ref. [50]). The c i denote the seven (complex) subtraction constants that are mandatory to obtain convergent dispersion integrals. These subtraction constants cannot be determined by dispersion theory alone and have to be fixed by more fundamental dynamical theories or, as performed here, by fits to experimental data. In particular, it might be possible to constrain their imaginary parts by three-body-unitarity considerations; this will be difficult in practice, however, as the threebody invariant mass is fixed to that of the decaying D-meson.
The constraint on the Mandelstam variables s + t + u = const. has allowed us to eliminate some of the subtraction constants, which are equivalent to the leading Taylor coefficients in an expansion around s/t/u = 0, for some of the single-variable functions in order to obtain a unique decomposition of the full decay amplitude. We have chosen the nonresonant I = 3/2 and I = 2 amplitudes for that purpose; see Ref. [26] for all details. Strictly speaking, the single-variable amplitudes in the above decomposition need an even higher degree of subtractions due to the high-energy behavior of the D-wave F 1/2 2 . The high-energy behavior of the decay amplitudes, and thus the single-variable amplitudes times their polynomial prefactors, is chosen to be consistent with the Froissart bound, which the F 1/2 2 amplitude cannot satisfy. 2 The thorough inclusion of the πK isospin 1/2 D-wave would thus necessitate more unknown subtraction constants that lower the predictability of the system. We therefore choose to include the F 1/2 2 amplitude heuristically in the sense that we exclude all crossed-channel projections of the F 1/2 2 amplitude in the representation above: the other (lower) partial waves are allowed to contribute to the inhomogeneityF 1/2 2 , but not vice versa. For an in-depth derivation of the full decay amplitude as well as the single-variable amplitudes we refer the reader to Ref. [26].
The set of dispersion relations (8) is linear in the singlevariable amplitudes F I J (s) as well as in the subtraction constants. To solve the system it is therefore beneficial to exploit this linearity and construct a basis of the solution space. The basis functions that span this solution space can be obtained by choosing a maximal set of linearly independent subtractionconstant configurations and solve the integral equations for each configuration. We choose the subtraction-constant configuration c j = δ i j , j ∈ {0, 1, . . . , n − 1} for the ith basis function M i (s, t, u). Thus the general solution can be written as a linear combination (and n = 7 in the system (8)). The explicit numerical solution strategy to determine the basis functions via matrix inversion is discussed in detail in Ref. [26].

Experimental comparison
In this section we perform fits of the dispersively determined decay amplitudes, displayed above, to the experimental 2 For the assumed asymptotic behavior of the input phase shifts that determines the one of the Omnès functions and hence the single-variable amplitudes, we refer to Ref. [26]. D + → K 0 S π 0 π + /K − π + π + data of the BESIII [39], CLEO [36], and FOCUS [35] collaborations. The D + → K 0 S π 0 π + Dalitz plot is totally dominated by the D + →K 0 π 0 π + decay, since the D + → K 0 π 0 π + decay channel is doubly Cabibbo-suppressed.
Previously, the same amplitudes have been compared to the D + → K − π + π + Dalitz plot data in detail [26]. Here we will focus on the D + → K 0 S π 0 π + data from the BESIII collaboration and subsequently perform combined fits to the D + → K 0 S π 0 π + /K − π + π + data sets to use the isospin relation between these channels to full capacity, as well as to check the extracted subtraction constants for consistency. Furthermore we will discuss the inclusion of the πK D-wave in more detail.

Comparison to the BESIII data
We begin with the comparison to the Dalitz plot data of the D + → K 0 S π 0 π + decay measured by the BESIII collaboration [39]. The experimental Dalitz plot, shown in Fig. 1 (upper plot), exhibits two prominent resonant structures, the ρ(770) and K * (892) resonances. Compared to the D + → K − π + π + Dalitz plot, which shows only a prominent K * (892) resonance, we have direct access to the ππ P-wave and therefore a much stronger constraint on the subtraction constants c 0 and c 1 . Since our dispersive approach requires elastic two-particle unitarity, Eq. (6), we restrict the fit region to pion-kaon two-particle energies below the η ′ K threshold, above which the major onset of inelasticities is seen phenomenologically [51][52][53], in particular in the S -wave. We can therefore assume elastic unitarity  to be a good approximation below this threshold. The provided data set comprises a binned t × u Dalitz plot with bin size 0.05 GeV 2 × 0.05 GeV 2 . We define the event distribution function of the efficiency-and background-corrected data by with (t i , u i ) being the center of the corresponding bin and 2δ = 0.05 GeV 2 the bin width. We perform two fit scenarios: without the πK D-wave F 1/2 2 (x) (fit 1) and including the D-wave (fit 2), analogously to Ref. [26]. The fit region is confined to t, s < (M η ′ + M K ) 2 and the χ 2 given by where N is the overall normalization of the corrected data with 746 bins in the considered fit region. Furthermore we define fit fractions as follows: where the P J (x) denote the angle-dependent prefactors of the corresponding single-variable amplitudes in the total amplitude. The fit results for the subtraction constants are summarized in Table 1, the emanating fit fractions are displayed in Table 2. In Fig. 1 (lower plot) we display the fitted theoretical Dalitz  (8). However, no region of particular disagreement is observed in the Dalitz plot. The moduli of the subtraction constants differ significantly from the ones extracted from D + → K − π + π + [26] and show smaller uncertainties. However, the phases of the subtraction constants attained in the BESIII fit are compatible (modulo 2π) with the CLEO/FOCUS phases of Ref. [26], with the exception of the phase of c 1 . We note that c 1 is the linear ππ P-wave subtraction constant, which contributes only indirectly via (chargeexchange) rescattering to the D + → K − π + π + decay amplitude. The phases of the F 1/2 0 subtraction constants mutually agree modulo π and can be chosen nearly real with an overall phase factored out, similar to what has been observed in Ref. [26]. However, this does not hold for the F 1 1 amplitude.
Strictly speaking, due to isospin symmetry, the CLEO, FO-CUS, and BESIII fits should result in the same values for the subtraction constants. Seeing that the BESIII fit clashes with the combined CLEO/FOCUS results of Ref. [26], it is however doubtful that a combined fit proves to be successful. With the overall χ 2 given by the sum of the individual χ 2 values, we attempt to perform a simultaneous fit of all three data sets (CLEO, FOCUS, and BESIII) available to us. The fit results in χ 2 combined /d.o.f. values of 1.7 ± 0.1 (2.5 ± 0.2) for fit 1 (fit 2): the inclusion of the πK D-wave in the combined fit considerably worsens the quality of the data description. This suggests that the heuristic inclusion of the πK D-wave, which is necessary to obtain sensible fit fractions in the CLEO fit [26], may seem sufficient for the individual fits, but is clearly not for a combined analysis.

Alternative D-wave model
In this section we want to assess the origin of the bad fit qualities in the combined analysis. Since the prime candidate is the only partially included πK I = 1/2 D-wave, we attempt to improve on the latter by oversubtracting the amplitude F 1/2 2 once, with the aim to obtain a more flexible description of the D-wave strength: . (14) Essentially this is equivalent to adding the πK D-wave Omnès function times a (complex) normalization constant c 7 to the original set of amplitudes. In principle, this worsens the highenergy behavior of the D-wave contribution even more; in practice, given our prescription to drop the D-wave contributions to the inhomogeneities, it does not render this inconsistency any more severe. We have previously justified the inclusion of the D-wave amplitude as in Eq.   the crossed-channel S -and P-waves, while neglecting the projections of F 1/2 2 onto the other partial waves-by alluding to low-energy processes such as those calculated in chiral perturbation theory [26]. For instance, the I = 0 ππ D-wave scattering length is almost completely given by crossed-channel dynamics [54,55], not by the low-energy tail of the f 2 (1270) resonance. Similarly, but outside the realm of applicability of chiral dynamics, the subleading F-wave in ω → 3π is dominated by crossed-channel amplitudes, not by the ρ 3 (1690) [22]. To assume that this picture could be extended beyond the nearthreshold region, up to energies where the D-wave becomes resonant, was clearly too optimistic. In this sense, the new subtraction constant c 7 could be linked to an independent coupling constant to the K 2 (1430) tensor resonance: it is the only parameter that allows to adjust the strength of the πK D-wave independently of the crossed-channel dynamics.
We again consider two fit scenarios: a fit to the BESIII data alone and a combined fit to the CLEO, FOCUS, and BESIII data sets. The results are summarized in Table 3, and the ensuing fit fractions in Table 4.
The BESIII fit results for the subtraction constants turn out to be very similar to the BESIII fits 1 and 2, see Table 1 for comparison. This is an anticipated result, since in the previous fit scenarios (fits 1 and 2) the inclusion of the D-wave did not change the subtraction constants substantially, but led to a poorer χ 2 result. Similarly the fit fractions, comparing the previous fit scenario 2 and the alternative D-wave BESIII fit, are alike, with the exception of the F 1/2 2 single-variable amplitude, which is slightly smaller. The fit quality with the more flexible D-wave is improved to χ 2 BES /d.o.f. = 1.08 ± 0.01 compared to χ 2 /d.o.f. = 1.27 ± 0.01 (χ 2 /d.o.f. = 1.35 ± 0.07) for fit 1 (fit 2) in the previous section.
With all data sets combined, we obtain a combined χ 2 /d.o.f. = 1.22 ± 0.03, which presents a huge improvement over the previous combined fits 1/2 (χ 2 /d.o.f. = 1.7 ± 0.1/2.5 ± 0.2). The χ 2 thus advocates that the discrepancy in fits 1/2 comes from the heuristically built-in D-wave. However, we note again that a thorough inclusion of the D-wave in the Khuri-Treiman formalism would necessitate further subtrac-tion constants in the other partial waves.
The subtraction-constant results are very similar to the individual BESIII fit values, but the χ 2 BES worsens from 1.08 ± 0.01 to 1.26 ± 0.04. In contrast, we find that the obtained χ 2 CLEO and χ 2 FOCUS values are similar to the individual fit qualities of Ref. [26], although the subtraction constants are significantly different. This suggests that the D + →K 0 π 0 π + data constrains the subtraction constants much better than the available D + → K − π + π + one. The fit fractions for the D + →K 0 π 0 π + decay amplitude are very similar to the BESIII fit while the large constructive interference effects seen in the CLEO/FOCUS fits of Ref. [26] do not show up as prominently in the D + → K − π + π + fit fractions. Additionally, we observe that the fit fractions of the non-resonant waves F 2 0 and F 3/2 0 reduce compared to the individual BESIII and CLEO/FOCUS fits.
Clearly, our χ 2 /d.o.f. = 1.22 ± 0.03 for the combined fit, though largely improved, is still not remarkably good, given the large number of data points. To put this into perspective, however, we wish to point out that most previous theoretical studies refrained from fitting Dalitz plots at all [27][28][29][30][31], while Ref. [38] performs a fit to pseudodata that does not allow for a sensible statistical interpretation. Isobar-model fits to the FO-CUS [35,37] and BESIII [39]  In Fig. 2, we show moduli and phases of the extracted single-variable amplitudes, comparing the results of the CLEO/FOCUS fits of Ref. [26] to fit 2 to the BESIII data as well as the combined fit to all three data sets including the alternative D-wave description. For the purpose of comparison, the overall normalizations of the total amplitudes are chosen such that in Fig. 2, the peak strengths of the K * (892) resonance in the I = 1/2 P-wave coincide in all fits. The main observation is that the new fits of the present study including the BESIII data, while consistent with the previous results, constrain most amplitudes much better; above all others, this rather obviously holds true for the ππ P-wave. An important result is that the D + → K S π 0 π + Dalitz plot confirms the conclusion of Ref. [26] concerning the phase motion of the I = 1/2 πK S -wave, whose phase rises much more quickly than the elastic scattering phase shift [56]. This has been observed before [34,37], but does not contradict Watson's final-state theorem [57]: the phase is modified due to rescattering with the third particle in the final state.
The πK D-wave merits a further comment, as the different fits in Fig. 2 display quite different magnitudes of F 1/2 2 : in particular the fit to the BESIII data alone yields a much larger contribution than the others. Remember that before the oversubtraction of Eq. (14), this is indirectly determined by the other partial waves. While the D + → K S π 0 π + data in the BESIII fit show rather little sensitivity to the D-wave, but strong constraints on the subtraction constants, in particular the fit fractions in the D + → K − π + π + data depend significantly thereon [26]. Only the more flexible form allows to reconcile both, with destructive interference between the D-wave generated through crossedchannel rescattering and the new subtraction constant c 7 .
Finally, the seemingly stark modifications of the phases of the small, nonresonant amplitudes must not be overstated: they occur in places where the moduli are close to zero, hence the full amplitudes are not severely affected.

Conclusion
In this letter we have resumed the study of strong finalstate interactions in the D + → K − π + π + /K 0 π 0 π + decays utilizing dispersion relations in the form of Khuri-Treiman-type equations. The resulting dispersion relations satisfy analyticity, unitarity, and crossing symmetry by construction and generate crossed-channel rescattering between the three final-state particles. We solely rely on ππ and πK scattering phase-shift input. The seven complex subtraction constants are fitted to the D + → K − π + π + /K 0 π 0 π + data from the CLEO, FOCUS, and BESIII collaborations. The requirement of elastic unitarity restricts the comparison to energies below the onset of major inelasticities. We have focused on two main fit scenarios: the individual fit to the D + →K 0 π 0 π + Dalitz plot from the BESIII collaboration, and a combined fit to all the data sets available to us. In both scenarios we have in particular studied the impact of the πK D-wave.
The individual fits to the BESIII D + →K 0 π 0 π + data show that the Dalitz plot in the elastic region can be described reasonably well with an overall χ 2 /d.o.f. of around 1.3. The subtraction constants are much more constrained than previously by D + → K − π + π + data alone, and much more stable with respect to inclusion of the πK D-wave, which does not improve the fit quality. In particular the direct sensitivity to the ππ Pwave has a large impact. We confirm once more that crossedchannel rescattering effects significantly shape the phases of the partial wave amplitudes; in particular the isospin 1/2 πK Swave phase shows a much steeper rise compared to the elastic scattering phase shift.
We find, however, that the values of the subtraction constants do not agree between the two differently charged final states, and a combined fit does not lead to satisfactory results. This tension could be traced back to the πK D-wave, which is not fully consistently included in our Khuri-Treiman system. Introducing an additional subtraction in the πK D-wave, hence allowing for an additional free parameter therein, we demonstrate that the description of the BESIII data improves to χ 2 /d.o.f. ≈ 1.1 in the individual fit. Consistency of all three data sets (and both charge channels) is reinstated, with an overall χ 2 /d.o.f. of about 1.2.
We have therefore demonstrated successful steps towards the application of Khuri-Treiman-type amplitude representations to three-body decays of charmed mesons, exploiting the consequences of isospin symmetry and the coupling of the various partial wave through rescattering in an optimal way. Further, systematic applications of this formalism should be investigated in the future. These should expand on the effects of inelasticities and coupled channels [21,58] in order to extend the dispersive description to the complete D-meson Dalitz plots, investigate the effects of three-body unitarity, and search for ways to cements  [26] in blue. The overall normalization is chosen such that the absolute values in the K * (892) peak in the I = 1/2 P-wave agree. Right column: Phase motion of the single-variable amplitudes (BESIII in red, alternative D-wave turquoise, CLEO/FOCUS in blue) and input scattering phases (black) in radiant. The phases are fixed to zero at the two-particle (ππ, πK) thresholds. Note that we obtain two separate solutions for the F 2 0 phase. The dashed lines visualize the fitted area; for the πK amplitudes from threshold to the η ′ K threshold and the full phase space for the ππ amplitudes. include higher partial waves consistently without a proliferation of free parameters.