Dalitz analysis of $D^{0}\to K^{-}\pi^{+}\eta$ decays at Belle

We present the results of the first Dalitz plot analysis of the decay $D^{0}\to K^{-}\pi^{+}\eta$. The analysis is performed on a data set corresponding to an integrated luminosity of 953 $\rm{fb}^{-1}$ collected by the Belle detector at the asymmetric-energy $e^{+}e^{-}$ KEKB collider. The Dalitz plot is well described by a combination of the six resonant decay channels $\bar{K}^{*}(892)^0\eta$, $K^{-}a_0(980)^+$, $K^{-}a_2(1320)^+$, $\bar{K}^{*}(1410)^0\eta$, $K^{*}(1680)^-\pi^{+}$ and $K_2^{*}(1980)^-\pi^{+}$, together with $K\pi$ and $K\eta$ S-wave components. The decays $K^{*}(1680)^{-}\to K^{-}\eta$ and $K_{2}^{*}(1980)^{-}\to K^{-}\eta$ are observed for the first time. We measure ratio of the branching fractions, $\frac{\mathcal{B}(D^{0}\to K^{-}\pi^{+}\eta)}{\mathcal{B}(D^{0}\to K^{-}\pi^{+})}=0.500\pm0.002{\rm(stat)}\pm0.020{\rm(syst)}\pm0.003{\rm (\mathcal{B}_{PDG})}$. Using the Dalitz fit result, the ratio $\frac{\mathcal{B}(K^{*}(1680)\to K\eta)}{\mathcal{B}(K^{*}(1680)\to K\pi)}$ is measured to be $0.11\pm0.02{\rm(stat)}^{+0.06}_{-0.04}{\rm(syst)}\pm0.04{\rm(\mathcal{B}_{\text{PDG}})}$; this is much lower than the theoretical expectations ($\approx1$) made under the assumption that $K^{*}(1680)$ is a pure $1^{3}D_1$ state. The product branching fraction $\mathcal{B}(D^0\to [K_2^{*}(1980)^-\to K^{-}\eta]\pi^{+})=(2.2^{+1.7}_{-1.9})\times10^{-4}$ is determined. In addition, the $\pi\eta^{\prime}$ contribution to the $a_0(980)^{\pm}$ resonance shape is confirmed with 10.1$\sigma$ statistical significance using the three-channel Flatt\'{e} model. We also measure $\mathcal{B}(D^0\to\bar{K}^{*}(892)^0\eta)=(1.41^{+0.13}_{-0.12})\%$. This is consistent with, and more precise than, the current world average $(1.02\pm0.30)\%$, deviates with a significance of more than $3\sigma$ from the theoretical predictions of (0.51-0.92)%.


I. INTRODUCTION
The understanding of hadronic charmed-meson decays is theoretically challenging due to the significant nonperturbative contributions, and input from experimental measurements thus plays an important role [1][2][3]. We present a Dalitz plot (DP) analysis [4] to study the dynamics of three body decay D 0 → K − π + η. This decay is Cabibbo-favored (CF) and proceeds via the c → sud transition. Because of isospin symmetry, intermediate states of this decay (e.g. excited kaon states decaying into Kπ or Kη), and a-family mesons decaying into πη, are similar to those in D 0 → K 0 S π 0 η. The DP analysis of the latter channel has previously been performed, and the intermediate channels K 0 S a 0 (980) 0 andK * (892) 0 η [5] were found to be dominant, but additional components of a non-resonant amplitude, K * 0 (1430)η, K 0 S a 2 (1320), κη, and combinations of these processes, were found to contribute significantly. However, the statistical power of that sample was too limited for precise measurements to be made. The D 0 →K * 0 η decay is sensitive to the W -exchange diagram, which is important for the theoretical understanding of charm decays. The theoretical predictions of the branching fraction of this mode vary in the range (0.51-0.92)% depending on the method [1][2][3]. This is consistent with, but smaller than, the current experimental result of (1.02 ± 0.30)% [6] obtained in the D 0 → K 0 S π 0 η final state [5]. A more precise measurement of this branching fraction from D 0 → K − π + η decays would test the theoretical predictions. The K * 0 (1430) ± → K ± η decay was observed by the BABAR experiment [7] and is awaiting confirmation. Experimentally, the ratio of K * 0 (1430) decaying into Kη and Kπ is 0.09 +0.03 −0.04 [6], which is consistent with the theoretical prediction of 0.05 [8] which was made with the assumption that it is a pure 1 3 P 0 state.
Decays of some other excited kaons to Kη, including K * (1410), K * (1680) and K * 2 (1980), were predicted by Refs. [8,9] but have not yet been observed. These states may have some interesting properties; K * (1410) may not be a simple 2 3 S 1 state, and K * (1410) and K * (1680) are predicted to be a mixture of the 2 3 S 1 and 1 3 D 1 states. Assuming K * (1410) and K * (1680) are pure 2 3 S 1 and 1 3 D 1 states, respectively, the relative branching ratio of K * (1680) to Kη and Kπ should be close to one (1.18 in Ref. [8] and 0.93 in Ref. [9]) Experimentally, no branching ratio measurement for the former channel has previously been made.
The nature of the a 0 (980) is still not clear. Since it is a dominant intermediate resonance in D 0 → K − π + η, we can collect a large sample of a 0 (980) + decays to study its character further, e.g. to confirm the πη ′ contribution to the a 0 (980) lineshape in a Flatté model as measured by BESIII [10]. Such a study can also help determine the πη and KK contributions to a 0 (980) precisely and understand its quark component.
Wrong-sign (WS) decays play an important role in studies of D 0 -D 0 mixing and CP violation such as the first observation of D 0 -D 0 mixing [11]. One possible mode for this, D 0 → K + π − η, will be reconstructed at Belle II, which aims at a data set fifty times [12] larger than that currently available from Belle. A time-dependent Dalitz analysis of this mode can be used to measure charm-mixing parameters, and for such a measurement an amplitude analysis of the right-sign decay, D 0 → K − π + η, is needed to obtain the CF decay model. This paper is organized as follows. Section II briefly describes the Belle detector and data samples, and Sec. III discusses event selection and parameterizations of signal and background and presents the measurement of the overall branching fraction. In Sec. IV, we report the results of the DP analysis. The evaluation of the systematic uncertainties are discussed in Sec. V. Further study and discussion of the Dalitz fit results are presented in Sec. VI. Finally, the conclusions are presented in Sec. VII.

II. BELLE DETECTOR AND DATA SETS
We perform a first Dalitz analysis of the decays D 0 → K − π + η [13] using 953 fb −1 of data collected at or near the Υ(nS) resonances (n=1, 2,3,4,5), where 74% of the sample is taken at the Υ(4S) peak, with the Belle detector [14] operating at the KEKB asymmetric-energy e + e − collider [15]. The Belle detector is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), a barrel-like arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter comprised of CsI(Tl) crystals located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-return located outside of the coil is instrumented to detect K 0 L mesons and to identify muons. A detailed description of the Belle detector can be found elsewhere [14].

III. EVENT SELECTION AND YIELDS
The signal decay chain consists of D * + → D 0 π + s with D 0 → K − π + η and η → γγ; D * mesons are produced in e + e − → cc processes, and the charge of the slow pion π s tags the flavor of the D 0 meson [16]. To ensure charged tracks are well reconstructed, each is required to have at least two associated hits of the SVD in the beam and azimuthal directions, separately. The slow-pion candidates are required to have the signed distances from the pivotal point to the helix to be within ±1.0 cm in the transverse plane and within ±3.0 cm along the direction opposite to the positron beam. A charged track is identified as a kaon by requiring a ratio of particle identification likelihoods [17] L K /(L K + L π ) > 0.7, which are constructed using CDC, TOF, and ACC information; otherwise the track is assumed to be a pion. This requirement has efficiencies of 85% and 98% and misidentification rates of 2% and 10% for kaons and pions, respectively. The photon candidates are reconstructed from ECL clusters unmatched to any charged track. The ratios of their energy deposits in a 3×3 array of CsI(Tl) crystals to that in a 5×5 array centered on the crystal with maximum deposited energy are required to be more than 0.8. The energies of photon candidates used to form η, E γ , must exceed 60 or 120 MeV in the barrel or endcap region. The η candidates must have γγ mass within −0.06

+0.05
GeV/c 2 of the nominal mass [6] which takes into account the asymmetric resolution, and to have momentum in the laboratory frame, p η , larger than 1 GeV/c. Furthermore, we require | cos θ η | = | Eγ1−Eγ2 Eγ1+Eγ2 | · Eη pη to be less than 0.8, which is optimized to suppress combinatorial background. A large set of simulated signal Monte Carlo (MC) samples, more than 25 million events, is produced to study the efficiency. These are generated uniformly in phase space with the Evtgen [18] and Jetset [19] software packages, and the detector response is modeled by the Geant3 [20]. The final-state radiation effect is taken into account using the Photos [21] package.
Kaon and pion tracks with opposite charge are required to form a common vertex (the D 0 decay position) with fit quality χ 2 v . The η candidates are mass-constrained assuming that they are produced at this decay vertex, and the resultant η momentum is added to the Kπ system to obtain the D 0 momentum. The invariant mass of Kπη, M , is required to satisfy the condition 1.80 GeV/c 2 < M < 1.92 GeV/c 2 . Then, the D 0 production vertex is constrained to the e + e − interaction point, with fit quality χ 2 b . The π s track is refit to this D 0 production vertex, with a fit quality denoted χ 2 s , to improve the resolution of the released energy in D * + decay, Q ≡ M Kπηπs − M Kπη − m πs . The value of Q is required to be less than 15 MeV/c 2 to suppress further combinatorial background. The D * momentum in the center-of-mass frame, p * (D * ), is required to be greater than 2.4, 2.5, or 3.1 GeV/c for data below, on, or above Υ(4S) energy, to reduce high-multiplicity events and combinatorial background. A consequence of this requirement is that the D 0 candidates from B decays are removed. After applying all of these selection criteria, there are on average 1.3 signal decay candidates per event. A best-candidate selection (BCS) method is applied to multi-candidate events, retaining as the best candidate the one with the smallest sum of vertex-fit qualities, A mass-constrained fit is then applied to the D 0 meson to improve resolution on the Dalitz variables, M 2 Kπ and M 2 πη . To extract yields of signal and background, a fit of the two-dimensional distribution of M and Q is performed. For the signal, the probability density function (PDF) in M is described by the sum of a double Gaussian and a double bifurcated Gaussian, with a common mean value (µ); the PDF in Q is described by the sum of a bifurcated Student function, a bifurcated Gaussian function, and a bifurcated Cruijff function [22], where the mean values and widths are correlated to the M value by a second-order polynomial function of |M − µ|. For a real signal D 0 combined with a random π s (named the random π s background), the M distribution uses the same PDF as for the signal and the Q distribution uses a threshold function, f (Q) = Q α e −βQ . This random background will be treated as signal, as it nearly consists of the same D 0 decay as the signal when the tiny fraction of DCS decay relative to CF decay is neglected. The combinatorial background is considered to have two components. A PDF smoothed by bilinear interpolation [23] is used for correlated combinatorial background, which has a correctly reconstructed π s from D * decay, but incorrectly reconstructed D 0 , whereas for other combinatorial background a third-order polynomial function of M and a threshold function of Q is used as a parameterization. The ratio between these two combinatorial backgrounds is fixed to that found using the generic MC. Figure 1 shows the M -Q combined fit for the experimental data. We obtain a signal yield of 105 197 ± 990 in the M and Q two-dimensional (2D) signal region of 1.85 GeV/c 2 < M < 1.88 GeV/c 2 and 5.35 MeV/c 2 < Q < 6.35 MeV/c 2 with a high purity (94.6 ± 0.9)%. These are the combinations that will be used for the fit to the Dalitz plot. To measure the branching fraction of the decay D 0 → K − π + η, we normalize the signal yield by the number of D 0 mesons produced in the decay D * + → D 0 π + . For normalization, we choose the D 0 → K − π + channel, which has a well-known rate of B = (3.950 ± 0.031)% [6]. We use the same selection criteria as are used D 0 → K − π + η but without the η. We extract the signal yield from the distribution of D 0 invariant mass in 1.78 GeV/c 2 < M < 1.94 GeV/c 2 and Q wide signal region |Q − 5.85| < 1.0 MeV/c 2 and find signal yields of 116 302 ± 510 for D 0 → K − π + η, and 2 597 343 ± 1 669 for D 0 → K − π + (with a high purity 98.3%) based on the Υ(4S) on-resonance data set. The efficiency ǫ(D 0 → K − π + η) = (5.34 ± 0.01)% and ǫ(D 0 → K − π + ) = (23.49 ± 0.02)% are determined based on Dalitz signal MC produced with the nominal Dalitz fit result shown in Table I for D 0 → K − π + η and signal MC for D 0 → K − π + . Taking into account the branching fraction B(η → γγ) = (39.41 ± 0.20)% [6], we find the ratio of branching fractions to be where the three uncertainties shown are statistical, systematic, and the uncertainty of branching fraction of η → γγ, respectively. Using the known D 0 → K − π + branching fraction, we measure the branching fraction where the last error is associated with uncertainty of the branching fractions of D 0 → K − π + and η → γγ. Many systematic uncertainties are canceled in the ratio measurement, and the dominant uncertainty is that of the η reconstruction efficiency (4%).

IV. DALITZ ANALYSIS
The isobar model [24] is applied for the amplitude of D 0 → (R → AB)C through a resonance R with spin-J (A, B and C are pseudoscalar particles). The decay amplitude is given by a coherent sum of individual contributions, consisting of a constant term a N R e iφNR for the non-resonant three-body decay, and different quasi-two-body resonant decays: Here m 2 AB and m 2 BC are Dalitz variables, and a R e iφR is a complex amplitude for the contribution of an individual intermediate resonance R. The amplitude and phase ofK * (892) 0 , having the largest fit fraction, are fixed to a K * (892) 0 = 1 and φ K * (892) 0 = 0. The matrix element M R for an intermediate resonant decay is given by where T R × Ω J is a resonance propagator. T R is a dynamical function for a resonance, described by a relativistic Breit-Wigner (RBW) with mass-dependent width, where p AB (p R ) is the momentum of either daughter in the AB (or R) rest frame, and M R and Γ R 0 are the nominal mass and width, Ω J describes the angular momentum that depends on the spin J by using the Zemach tensor [25,26], and F D and F R are Blatt-Weisskopf centrifugal barrier factors [27,28], describing the quark structure of the D 0 meson and intermediate resonance. The parameter of meson radius, R, is set to 5.0 (GeV/c) −1 and 1.5 (GeV/c) −1 for the D 0 meson and the intermediate resonances, respectively [26]. For the a 0 (980) contribution description, we use the Flatté formalism with three coupled channels, πη,K 0 K and πη ′ [10] where √ s is the invariant mass of πη; g i and ρ i are coupling constants and phase-space factors, respectively. For . The generalized LASS model [29,30] is used to parameterize the Kπ and Kη S-wave contributions: where s is the invariant mass squared of the Kπ or Kη system, q is the momentum of K in the Kπ or Kη rest frame, and δ B and δ R are phase angles of the non-resonant component and K * 0 (1430) component, respectively. They are defined as tan(δ R ) = M r Γ(m ab )/(M 2 r − m 2 ab ) and cot(δ B ) = 1/(aq) + rq/2, where a, r, B, φ B and φ R are real parameters and may be determined by amplitude analysis.
The DP fit is performed by an unbinned maximum likelihood method with where n is the number of D 0 candidates in the M and Q 2D signal region and f i s is the event-by-event fraction of signal obtained from the M -Q fit; the combinatorial background function, P b , is a smoothed PDF [23], determined from the DP in the M sideband region (1.755 GeV/c 2 < M < 1.775 GeV/c 2 or 1.935 GeV/c 2 < M < 1.955 GeV/c 2 ) and the Q signal region (5.35 MeV/c 2 < Q < 6.35 MeV/c 2 ). The signal PDF, P s , is calculated taking the reconstructionefficiency dependence on the Dalitz-plot variables into account, and normalized in the Dalitz plot region.
This efficiency distribution ǫ(m 2 Kπ,i , m 2 πη,i ) is obtained from a high-statistic signal-MC sample and takes into account the known difference in particle identification efficiency for charged tracks between MC and data. These correction factors depend on the momentum and polar angle of individual charged track. The fit fractions (FF) of each intermediate component are calculated across the DP region as The FF uncertainties are evaluated using a Toy MC method in which the sampling takes into account the considerations between the amplitudes and phases by propagating the full covariance matrix obtained by the DP fit.
Fifteen possible intermediate resonances [31] were initially considered in the Dalitz analysis. We foundK * (1410) 0 andK * (1680) 0 have a phase-angle difference of approximately 180 • and similar behavior in the DP, therefore, it is hard to separate them. In order to ensure stability of the fit to the Dalitz plot, onlyK * (1410) 0 is kept, while a possibleK * (1680) 0 contribution is considered as a source of systematic uncertainty. Therefore for the rest of this paper, K * (1410) 0 represents the contribution of K * (1410) 0 , K * (1680) 0 Fig. 2 (b-d). The statistical significance of each component is larger than 10σ. In particular, the statistical significance of the Kη S-wave component with K * 0 (1430) − is greater than 30σ, and K * (1680) − → K − η and K * 2 (1980) − → K − η are observed for the first time and have statistical significances of 16σ and 17σ, respectively. The fitted magnitudes and phases of intermediate components are listed in Tab. I, together with corresponding fit fractions, where statistical uncertainties are obtained from 500 sets of toy MC samples, and systematic uncertainties take into account model uncertainties and other systematic uncertainties as discussed in Sec. V. The fact that the sum of fit fractions is greater than 100% indicates significant destructive interference. Table II shows the fitted parameters of LASS model in Eq. (7) for the Kπ and Kη S-wave components. Various Dalitz models, including the nominal model used in the fit to final experimental data, are produced using MC to perform tests for any possible bias, and to check that the input and output Dalitz parameters are consistent. We also checked for the existence of possible multiple solutions in the fit, with likelihood scanning of each of the free parameters. In addition, 100 sets of Dalitz fits were performed by sampling the initial values of free parameters uniformly in an interval around their final values. No multiple solutions were found.  To investigate the parameters of the Flatté formulation of the a 0 (980) + lineshape, the Dalitz fit based on the nominal model with free g π ± η is also performed and this yields g π ± η = 0.596 ± 0.008(stat) GeV/c 2 . This value is consistent with the measurement of BESIII, 0.607 ± 0.011 GeV/c 2 [10]. The significance of the πη ′ contribution is tested and the results with floated g πη ′ and fixed g πη ′ = 0 give ∆(−2 ln L) = 102 with ∆(d.o.f) = 1, which indicates a πη ′ contribution with 10.1σ statistical significance. The fitted g πη ′ = 0.408 ± 0.018(stat) GeV/c 2 is also consistent with the BESIII measurement of g πη ′ = 0.424 ± 0.050 GeV/c 2 [10].   [6]. To account for the Kπ and Kη S-wave components, the model used in the fit is modified by adding a wide resonance κ described by a complex pole function [32] for a Kπ S-wave, and K * 0 (1950) − described by RBW for a Kη S-wave. The nonsignificant resonance a 0 (1450) + is added to evaluate the πη S-wave component uncertainty. We also use aK * (1680) 0 resonance instead of aK * (1410) 0 contribution.
The systematic uncertainty due to the Dalitz distribution of combinatorial background is evaluated by (1) varying the M sideband region within a shift of ±5 MeV/c 2 , and by (2) correcting the Dalitz distribution of experimental data in the M sideband by the ratio of combinatorial background in the M signal and sideband regions from generic MC. The larger difference is assigned as the systematic uncertainty due to the background distribution. The systematic uncertainty related to efficiency is estimated in two ways: (1) removing the correction for PID efficiency, and (2) shifting the p * (D * ) limit by ±0.05 GeV/c to consider possible discrepancy between MC and experimental data in p * (D * ) spectrum. These uncertainties are combined quadratically to give a systematic uncertainty due to efficiency. Comparing with the nominal fit model, the difference in the fit results when the signal fraction is varied by ±1σ (as determined from the M -Q fit) is taken as the systematic uncertainty due to the uncertainty in the fraction of signal in the sample. A shift of the signal region by ±5 MeV/c 2 in M or ±0.1 MeV/c 2 in Q is applied to estimate the effect of the signal region selection. The larger difference in fit fraction is kept as the uncertainty due to this source. The uncertainty of multi-candidate selection is estimated by randomly selecting one of the multi-candidates as the best candidate instead of our nominal BCS method. The sources of systematic uncertainty considered are summarized in Tab. III. Individual uncertainties are added in quadrature.