Study of the decay $D^+\to K^*(892)^+ K_S^0$ in $D^+\to K^+ K_S^0 \pi^0$

Based on an $e^{+}e^{-}$ collision data sample corresponding to an integrated luminosity of 2.93 $\mathrm{fb}^{-1}$ collected with the BESIII detector at $\sqrt{s}=3.773 \mathrm{GeV}$, the first amplitude analysis of the singly Cabibbo-suppressed decay $D^{+}\to K^+ K_S^0 \pi^0$ is performed. From the amplitude analysis, the $K^*(892)^+ K_S^0$ component is found to be dominant with a fraction of $(57.1\pm2.6\pm4.2)\%$, where the first uncertainty is statistical and the second systematic. In combination with the absolute branching fraction $\mathcal{B}(D^+\to K^+ K_S^0 \pi^0)$ measured by BESIII, we obtain $\mathcal{B}(D^+\to K^*(892)^+ K_S^0)=(8.69\pm0.40\pm0.64\pm0.51)\times10^{-3}$, where the third uncertainty is due to the branching fraction $\mathcal{B}(D^+\to K^+ K_S^0 \pi^0)$. The precision of this result is significantly improved compared to the previous measurement.

Based on an e + e − collision data sample corresponding to an integrated luminosity of 2.93 fb −1 collected with the BESIII detector at √ s = 3.773 GeV, the first amplitude analysis of the singly Cabibbo-suppressed decay D + → K + K 0 S π 0 is performed. From the amplitude analysis, the K * (892) + K 0 S component is found to be dominant with a fraction of (57.1 ± 2.6 ± 4.2)%, where the first uncertainty is statistical and the second systematic. In combination with the absolute branching fraction B(D + → K + K 0 S π 0 ) measured by BESIII, we obtain B(D + → K * (892) + K 0 S ) = (8.69 ± 0.40 ± 0.64 ± 0.51) × 10 −3 , where the third uncertainty is due to the branching fraction B(D + → K + K 0 S π 0 ). The precision of this result is significantly improved compared to the previous measurement. This result also differs from most of theoretical predictions by about 4 σ, which may help to improve the understanding of the dynamics behind.

I. INTRODUCTION
The study of CP violation (CPV) in hadron decays is important for the understanding of the matterantimatter asymmetry in the universe.
In the charmed-meson sector, CPV effects in singly Cabibbosuppressed (SCS) D-meson decays are usually much larger than those in Cabibbo-favored and doubly Cabibbosuppressed decays. In 2019, the LHCb collaboration first observed a CPV effect in a combined analysis of D 0 → π + π − and D 0 → K + K − decays [1]. However, theoretical predictions of this CPV effect suffer from large variations compared to those in the K or B meson systems, mainly due to the large uncertainty in describing the non-perturbative dynamics in QCD in the charm region [2,3].
In recent years, the branching fractions (BFs) and CPV in the two-body hadronic decays of D → P P and D → V P have been studied in different QCDderived models [4][5][6][7], where P and V denote pseudoscalar and vector mesons, respectively. Generally, these theoretical calculations are in good agreement with experimental results, except for those of the SCS decay D + → K * (892) + K 0 S , whose amplitude consists of color-favored tree diagrams, W -annihilation diagrams, and penguin diagrams [7]. The topological diagrams can be found in Fig. 1. The measured and predicted values of B(D + → K * (892) + K 0 S ) are listed in Table I. The E687 collaboration reported a BF ratio of B(D + →K * (892) + (→K + π 0 )K 0 S ) B(D + →K 0 S π + ) = 1.1 ± 0.3 ± 0.4 [8]. This results in B(D + → K * (892) + K 0 S ) = (17 ± 8) × 10 −3 when combined with the world average of B(D + → K 0 S π + ) [9]. Although the predicted values are consistent with the experimental results, the experimental precision needs to be improved. A precise measurement of B(D + → K * (892) + K 0 S ) will provide a more stringent test of the theoretical models and help to deepen our understanding of the dynamics of charmed meson decays. Especially, this will enhance the predictive power on the CPV in charmed meson decays.
In this paper, the first amplitude analysis of the SCS decay D + → K + K 0 S π 0 is reported, based on a sample of D + D − pairs from e + e − collisions at a center-of-mass energy of √ s = 3.773 GeV corresponding to an integrated luminosity of 2.93 fb −1 [10,11] collected with the BESIII detector [12] at the Beijing Electron Positron Collider (BEPCII) [13]. At the energy √ s = 3.773 GeV, the pairs of D + D − are produced near threshold without any accompanying hadron [14]. Previously, the BESIII collaboration has measured the BF of D + → K + K 0 S π 0 to be (5.07 ± 0.19 ± 0.23) × 10 −3 [15]. In combination with the amplitude analysis results presented in this paper, the BF of D + → K * (892) + K 0 S can be determined with much improved precision. Charge conjugation is implied throughout the text.  [4], the factorization-assisted topologicalamplitude (FAT) approach with ρ-ω mixing [5], the topological diagram approach with only tree level amplitude (denoted as TDA[tree]) [6], and including QCD-penguin amplitudes (denoted as TDA[QCD-penguin]) [7]. For comparison, the previous experimental result [8,9] is also listed.

II. BESIII EXPERIMENT AND MONTE CARLO SIMULATION
The BESIII detector records symmetric e + e − collisions provided by the BEPCII storage ring, which operates with a peak luminosity of 1 × 10 33 cm −2 s −1 in the centerof-mass energy range from 2.0 to 4.9 GeV. BESIII has collected large data samples in this energy region [2]. The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identification modules interleaved with steel. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the ionization energy loss dE/dx resolution is 6% for electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps. More detailed descriptions can be found in Refs. [12,13].
Simulated data samples produced with a Geant4based [16] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector [17,18] and the detector response, are used to estimate background contributions and obtain the reconstruction efficiency. The simulation models the beam energy spread and initial state radiation (ISR) in the e + e − annihilations with the generator kkmc [19]. The 'inclusive MC sample' includes the production of DD pairs (including quantum coherence for the neutral D channels), non-DD decays of the ψ(3770), ISR production of the J/ψ and ψ(3686) states, and continuum processes which are incorporated in kkmc [19]. Known decay modes are modeled with evtgen [20] using the BFs published by the Particle Data Group (PDG) [9], and the remaining unknown charmonium decays are modeled with lundcharm [21]. The final state radiation from charged final state particles is incorporated using photos [22]. The inclusive MC sample is used to study background contributions and to estimate signal purity. In this work, two sets of signal MC samples are used. One sample is generated with a uniform distribution in phase space (PHSP) for the decay D + → K + K 0 S π 0 , called the 'PHSP MC sample', which is used to extract the detection efficiency maps along the Dalitz plot coordinates. The other sample is generated based on the fitted amplitudes from the amplitude analysis, called the 'DIY MC sample'. It is used to evaluate the fit quality and estimate the systematic uncertainty. The recoiling D − in these two sets of MC samples is forced to decay into six tag modes, discussed in Sec. III A.

III. EVENT SELECTION
Taking advantage of the threshold production of the D + D − sample, this analysis uses a double-tag method, which is illustrated in the following.
Charged particle tracks are reconstructed using the information of the MDC, and are required to have a polar angle θ with respect to the z-axis, defined as the symmetry axis of the MDC, satisfying |cos θ| < 0.93 and to have a distance of closest approach to the interaction c s d d point (IP) smaller than 10 cm along the z-axis (V z ) and smaller than 1 cm in the perpendicular plane (V r ). Those tracks used in reconstructing K 0 S → π + π − decays are exempt from these selection criteria. Particle identification (PID) for charged particle tracks is implemented by using combined information from the flight time measured in the TOF and the dE/dx measured in the MDC to form a PID probability L(h) for each hadron (h) hypothesis with h = π, K. Charged tracks are identified as pions when they satisfy L(π) > L(K), and as kaons otherwise.
Photon candidates from π 0 decays are reconstructed from the electromagnetic showers detected in the EMC crystals. The deposited energy is required to be larger than 25 MeV in the barrel region with |cos θ| < 0.80 and larger than 50 MeV in the end cap region with 0.86 < |cos θ| < 0.92. To further suppress fake photon candidates due to electronic noise or beam background, the measured EMC time is required to be within 700 ns from the event start time. To reconstruct π 0 candidates, the invariant mass of a photon pair is required to satisfy 0.115 < M γγ < 0.150 GeV/c 2 . To further improve the momentum resolution, the invariant mass of the photon pair is constrained to the nominal π 0 mass [9] by applying a one-constraint kinematic fit. The updated momentum of the π 0 is used in the further analysis.
K 0 S candidates are reconstructed through the decay K 0 S → π + π − by combining all pairs of oppositely charged tracks, without applying the PID requirement. These tracks need to satisfy |cos θ| < 0.93 and V z < 20 cm while no V r requirement is applied. A vertex fit is applied to pairs of charged tracks constraining them to originate from a common decay vertex, and the χ 2 of this vertex fit is required to be less than 100. The invariant mass of the π + π − pair needs to satisfy 0.487 < M π + π − < 0.511 GeV/c 2 . Here, M π + π − is calculated with the pions constrained to originate at the decay vertex.
To identify the tagged candidates, two kinematic variables, the beam-constrained mass M tag BC , and the energy difference ∆E tag , are defined as and where p tag and E tag are the three-momentum and energy of the tagged D − candidate in the rest frame of the initial e + e − collision system, and E beam is the beam energy. For a correctly reconstructed tagged candidate, M tag BC and ∆E tag are expected to be consistent with the nominal D − mass [9] and zero, respectively. In this work, the same ∆E tag requirements as those in a previous BESIII analysis [23] are used, which are listed in Table II. In each event, only the combination with the smallest |∆E tag | is kept for each tag mode. The tagged candidates are required to be within the region 1.863 < M tag BC < 1.879 GeV/c 2 .

B. Signal candidate selection
Signal candidates for D + → K + K 0 S π 0 are formed using the remaining tracks recoiling against the tagged D − .
Besides the selection requirements for charged and neutral tracks described in Sec. III A, some additional criteria are applied to improve the signal-to-background ratio.
For the K 0 S candidates, an additional secondary vertex fit is applied where the momentum of the reconstructed K 0 S candidate is constrained to be aligned with the direction from the IP to the K 0 S decay vertex, and the resulting flight length L is required to be larger than twice its uncertainty σ L . The χ 2 of the secondary vertex fit is required to be less than 500. For the π 0 candidates, the χ 2 of the kinematic fit is required to be less than 20.
To further identify the signal D + candidates, the energy difference, ∆E sig , is defined as where E sig is the energy of the signal D + candidate in the rest frame of the initial e + e − collision system. In each event, only the combination with the least |∆E sig | is kept as a signal candidate. The ∆E sig distributions of data and inclusive MC sample are shown in Fig. 2(a). We require −0.03 < ∆E sig < 0.02 GeV for signal candidates to be kept. To further improve the momentum resolution of the signal final state K + K 0 S π 0 , an additional kinematic fit is applied constraining the invariant mass of the signal final state to the nominal D + mass [9], and the total four-momentum of all reconstructed particles to the initial e + e − collision four-momentum. The updated four momenta are used for further analysis.
The recoil mass M rec is defined as where p e + e − is the e + e − collision initial four-momentum and p D + is the four-momentum of the D + signal candidate. The distribution of M rec in data and in the inclusive MC sample is shown in Fig. 2(b). The candidate events within 1.865 < M rec < 1.877 GeV/c 2 are selected and the signal purity is determined to be (97.4 ± 0.2)% according to the inclusive MC sample. After imposing all above selection criteria, the number of the signal events in data is measured to be 692.

IV. AMPLITUDE ANALYSIS
Only spin-zero particles are involved in the signal process D + → K + K 0 S π 0 , thus only two degrees of freedom are needed to describe the full kinematics. In the amplitude analysis, we choose the two Dalitz plot variables

A. Isobar model
The full amplitude of the decay process is described by the Isobar model [24], which is given by where the term c i = a i e iφi consists of the magnitude a i and the phase φ i for the specific intermediate process i. The amplitude A i denotes the decay amplitude of the process i, which is modeled in a quasi-two-body decay D → Cr, r → AB. Here, r is the possible resonance decaying into AB, and A, B and C each denote one of the final state particles K + K 0 S π 0 . A i is formulated as where W i is the spin factor, F r and F D are the Blatt-Weisskopf barrier factors [25], and T i is the dynamical function describing the intermediate resonance, as illustrated below.

Spin factor
The spin factor W i for the process D → Cr, r → AB is constructed based on the Zemach formalism [26]. The amplitudes for resonances with angular momenta larger than two are not considered due to the limited phasespace. The spin factor is expressed as where a 1 , a 2 , and a 3 are given by Here

Blatt-Weisskopf barrier factors
The Blatt-Weisskopf barrier factors F D and F r attempt to model the underlying quark structure of the parent particle in the decay D → rC and the subsequent decay r → AB, respectively. The expressions for the barrier factors, which are shown in Table III, are taken from Ref. [25]. In these expressions, p is the decay momentum of the particle A (C) in the rest frame of the mother particle r (D), while q is the decay momentum of the particle A (C) in the rest frame of the mother particle r (D) when the resonance is fixed at the corresponding nominal mass. The radii of D + and the intermediate resonance r are chosen as R D = 5 GeV −1 and R r = 1.5 GeV −1 , respectively [24,27].

Dynamical function
The dynamical function describes the line-shape of the intermediate resonance. For the specific process r → AB,  [25] under different angular momenta L.
the dynamical function is chosen as a relativistic Breit-Wigner function for most of the resonances and is written as where M r is the nominal resonance mass, M AB is the invariant mass of AB and Γ(M AB ) is the mass-dependent width defined as where Γ r is the nominal resonance width, J is the spin of the resonance and p AB and p r are the breakup momenta at M AB and M r , respectively. For the dynamical function of the Kπ S-wave, we choose the LASS parametrization [28]. It includes both the K * 0 (1430) resonances and a non-resonant part. We denote the full LASS parametrization as (K + π 0 ) S−wave or (K 0 S π 0 ) S−wave in the following. The LASS parametrization can be expressed as where F NR Kπ (φ NR Kπ ) and F Kπ ) are the magnitudes (phases) for the non-resonant and K * 0 (1430) components, respectively, and p Kπ is the breakup momentum of the Kπ system. The phase-shifts δ NR Kπ and δ K * 0 Kπ are defined as and where a scat is the scattering length, r eff is the effective interaction range, and M K * 0 is the nominal mass of the K * 0 (1430). The LASS parametrization corresponds to a K-matrix approach [29] describing a rapid phase shift coming from the resonant term and a slowly rising shift governed by the non-resonant term, with relative strengths F K * 0 Kπ and F NR Kπ . In the nominal fit, the LASS parameters are fixed according to the values measured by the BaBar and Belle collaborations [30], which are listed in Table IV.  [30].

B. Likelihood function
In the amplitude analysis, a maximum likelihood fit is performed by minimizing the negative log-likelihood (NLL), which is constructed on the Dalitz plot plane as where (x, y) denote the Dalitz plot coordinates y) is the efficiency function based on the smoothed histogram (following the method in Ref. [31]) from the PHSP MC sample, where the Dalitz plot and the corresponding projections of the PHSP MC are illustrated in Fig. 3, A i (x, y) is the decay amplitude of the i-th component in Eq. (6), c i is the free complex coefficient of the i-th component, and I ij is the normalization integral, which is defined as Here, the integral is calculated numerically by dividing the Dalitz plot plane into a grid of 3500 × 3500 square cells. No background contribution is included in the NLL in the nominal fit, exploiting the high signal purity.

C. Fit fraction and goodness-of-fit test
The fit fraction (FF) f i for the i-th component is calculated with where the integral is calculated using the same numerical method as for the integral in Eq. (15). Note that the sum of the FFs is not necessarily equal to unity due to the interferences between different components. To obtain the corresponding statistical uncertainties, the values of the fitted coefficients c i are randomly modified for 1000 times according to the information of the covariance matrix and the root-mean-square values of the distributions of the modified FFs are taken as the statistical uncertainties.
To examine the quality of the nominal fit, goodness-offit tests on three different projections of the Dalitz plot are performed using the fitted results. When calculating χ 2 for each projection, an adaptive binning is adopted to ensure that the minimum number of events is larger than 10 to obey the Gaussian assumption. The χ 2 value is calculated by using the number of events in data n k data and the expected number n k fitted from the nominal fit in the k-th bin and is formulated as

V. FIT RESULTS
To perform the amplitude analysis, the open-source framework GooFit [32] is used accelerating the fit speed using parallel processing computing. In the fit, the resonance K * (892) + is chosen as the reference, whose phase and magnitude are fixed to be one and zero, respectively. First, the fit of data is performed with the amplitudes containing K * (892) + andK * (892) 0 , which are clearly observed in the corresponding invariant mass spectra. Then, two S-wave Kπ components, (K + π 0 ) S−wave and (K 0 S π 0 ) S−wave are included. The statistical significances of these two components, calculated by the change of the log-likelihood values ∆(NLL) with and without including the component and taking into account the change of the number of degrees of freedom, are both found to be larger than 5 σ. Besides these four components, in addition K * (1410), K * 2 (1430), a 0 (980), a 0 (1450), ρ(1450), ρ(1700) and (Kπ) P−wave components were also tested, but their statistical significances are all lower than 5 σ, thus they are not included in the nominal fit. Here, the (Kπ) P−wave denotes the non-resonant P−wave contribution, which is modeled with the same as vector resonances like K * (892) + while the dynamical function is set to be constant.

VI. SYSTEMATIC UNCERTAINTY
The systematic uncertainties for the resonance amplitudes D + → K * (892) + K 0 S and D + →K * (892) 0 K + are discussed below. While the two S-wave Kπ components are included in our fit in order to improve the fit quality, in general the limited statistics of the data sample does not allow for detailed studies on these contributions so that we limit our systematic studies to the D + → K * (892) + K 0 S and D + →K * (892) 0 K + results. They are categorized into the following sources: (I) amplitude components, (II) input parameters for resonances, (III) radius of the meson (R r and R D ), (IV) background, (V) fit bias and (VI) efficiency. The results of the systematic uncertainties for phases and FFs are summarized in Table VII,   Source nominal fit with randomly added additional components from the list K * (1410), K * 2 (1430), a 0 (980), a 0 (1450), ρ(1450), ρ(1700), and (Kπ) P−wave . The magnitude of the additional component is randomly distributed in the range between zero and the fitted magnitude of the D + →K * (892) 0 K + component, and its phase is randomly distributed in the range [0, 2π). Then, the fit procedure is repeated for each toy MC sample. From these fits, we obtain pull distributions for the fit results φK * (892) 0 , FF K * (892) + , and FFK * (892) 0 compared to the nominal result that are well described by a Gaussian. The widths of the corresponding Gaussian functions describing the pull distributions are assigned as the associated systematic uncertainties.
(II) Input parameters: In the nominal fit, the masses and widths of K * (892) + andK * (892) 0 are fixed to the values in the PDG [9] and the parameters of the LASS model are fixed according to Ref. [30]. To estimate the corresponding systematic uncertainties, the fit procedure is repeated by varying each of the fixed parameters by ±1σ. The quadratic sum of the maximum relative variations for each parameter is taken as the systematic uncertainty.
(III) Radii of the mesons: To estimate the relevant systematic uncertainties originating from fixing the R D value, the fits are performed with alternative R D values between 3 GeV −1 and 7 GeV −1 . The maximum relative variations of the fit results are taken as the relevant systematic uncertainties. In case of the R r value, which is one of the dominant sources of systematic uncertainty, toy MC samples are generated based on the fit results with randomly distributed R r parameters in the range[0, 3] GeV −1 . These toy MC samples are then fitted with the R r parameter fixed to the default val-ue of 1.5 GeV −1 . The observed pull distribution can be described by a Gaussian function. The width of the Gaussian function used to fit the pull distribution is taken as systematic uncertainty. The quadratic sum of these two uncertainties is taken as the corresponding systematic uncertainty.
(IV) Background: In the nominal fit, the background is neglected due to high signal purity. To estimate the associated systematic uncertainty, the NLL is alternatively constructed as where f is the signal fraction calculated from the M rec distribution in the inclusive MC sample and B(x, y) is the background distribution modeled with a smoothed histogram [31] constructed from the inclusive MC sample. After minimizing the NLL in Eq. (18) with the components of the nominal solution, the relative variation of the fit results is found to be smaller than 1% of the corresponding statistical uncertainty. The variations are assigned as the systematic uncertainties.
(V) Fit bias: To understand the potential effect of the fit bias, a series of DIY MC samples with same statistics as in data are generated. The fit procedure is repeated for each DIY MC sample and the pull distributions compared to the nominal fit result is obtained. Any deviation from a mean of zero of the pull distribution is assigned as a systematic uncertainty.
(VI) Efficiency: Uncertainties from the efficiencies for charged particle tracking and PID, as well as the reconstruction of K 0 S and π 0 candidates have been studied with different control samples in previous works, see Refs. [33][34][35] for examples. To estimate the corresponding systematic uncertainties, we use correction factors comparing the estimated efficiencies in data and MC simulation, ε Data /ε MC , to re-weight the efficiency function. We obtain a modified efficiency η ′ (x, y), which is then used instead of η(x, y) in Eqs. (14) and (15). The relative variation of the fit results using the modified efficiency is taken as the corresponding systematic uncertainty. Additionally, the systematic uncertainties due to the M rec and ∆E requirements are studied by slightly shifting the boundaries within 1 MeV/c 2 and 1 MeV, respectively. The modified efficiency functions are obtained and the fit procedure is repeated. Finally, the quadratic sum of the above variations is assigned as the systematic uncertainty related to the efficiency.

VII. SUMMARY
To summarize, based on an e + e − collision sample corresponding to an integrated luminosity of 2.93 fb −1 collected with the BESIII detector at √ s = 3.773 GeV, the first amplitude analysis of D + → K + K 0 S π 0 is carried out.