Study of the Decay $D^{+}_{s}\rightarrow \pi^{+}\pi^{+}\pi^{-}\eta$ and Observation of the W-annihilation Decay $D^{+}_{s}\rightarrow a_0(980)^+\rho^0$

The decay $D^{+}_{s}\rightarrow \pi^{+}\pi^{+}\pi^{-}\eta$ is observed for the first time, using $e^+e^-$ collision data corresponding to an integrated luminosity of 6.32 fb$^{-1}$, collected by the BESIII detector at center-of-mass energies between 4.178 and 4.226 GeV. The absolute branching fraction for this decay is measured to be $\mathcal{B}(D^+_s \to \pi^+ \pi^+ \pi^- \eta) = (3.12\pm0.13_{\rm stat.}\pm0.09_{\rm syst.})$%. The first amplitude analysis of this decay reveals the sub-structures in $D^{+}_{s}\rightarrow \pi^{+}\pi^{+}\pi^{-}\eta$ and determines the relative fractions and the phases among these sub-structures. The dominant intermediate process is $D^{+}_{s}\rightarrow a_1(1260)^+ \eta, a_1(1260)^+ \rightarrow \rho(770)^0\pi^+$ with a branching fraction of $(1.73 \pm 0.14_{\rm stat.} \pm 0.08_{\rm syst.})$%. We also observe the W-annihilation process $D^{+}_{s}\rightarrow a_0(980)^+\rho(770)^0$, $a_0(980)^+ \to \pi^+ \eta$ with a branching fraction of $(0.21\pm0.08_{\rm stat.}\pm0.05_{\rm syst.})$%, which is larger than the branching fractions of other measured pure W-annihilation decays by one order of magnitude.


PACS numbers:
Since the discovery of the charm quark in 1974, the hadronic decays of charmed mesons have been extensively studied both experimentally and theoretically.However, making precise Standard Model predictions for exclusive weak charm meson decays is rather difficult, because the charm quark mass is not heavy enough to allow for a sensible heavy quark expansion, as corrections of higher orders in 1/m c become very important, nor is it light enough for the application of chiral perturbation theory [1].Studies of hadronic D + s decays provide insight into its decay mechanisms.The Cabibbo-favored (CF) hadronic D + s decays mediated via a c → sW, W → u d transition, producing states with hidden strangeness (ss quarks in η), have relatively large branching fractions (BFs), although some of them have still not been measured.The Particle Data Group (PDG) [2] reports that the missing hadronic decays of D + s with η in the final state contribute a fraction of (7.1 ± 3.2)%.Among them, the D + s → π + π + π − η decay is expected to have a large BF, but until now it has not been observed.
The topological diagram analysis of the hadronic D + s decays [3] reveals the importance of weak annihilation contributions [4].The amplitude of the long-distance weak annihilation induced by final-state interactions may be comparable to the tree amplitude [1,4].The BF of the pure W-annihilation (WA) process with a P P final state, D + s → π 0 π + , is less than 0.037% [5], and that of the V P final state, D + s → ρ(770) 0 π + , is 0.019% [2], while the BF of the WA process with an SP final state, D + s → a 0 (980) +(0) π 0(+) , a 0 (980) +(0) → π +(0) η, is about 1.46% [7].Here, V , P , and S denote vector, pseudoscalar, and scalar mesons, respectively.Since the direct production of the a 0 (980)π system via the cs → W + → u d → a 0 (980)π transition violates G-parity conservation [8], the WA process D + s → a +(0) 0 π 0(+) is suppressed.The origin of the abnormally large BF for the SP decay mode is still controversial.Reference [9] argues that the SP decay mode can be produced via internal emission due to which a 0 (980) can be dynamically generated from K K finalstate interactions in coupled channels.Reference [10] claims that D + s → a 0 (980) +(0) π 0(+) receives the main contribution from D + s → ρ(770) +(0) η through a triangle rescattering process.To date, there have been no theoretical or experimental studies of the WA process with a V S final state.The process D + s → a 0 (980) + ρ(770) 0 can proceed via a WA transition with G-parity conservation.The decay diagram for this decay is shown in Fig 1 .A measurement of this process is critical to distinguish various weak annihilation mechanisms.In addition, the underlying structure of the resonance a 0 (980) + plays an important role in the decay mechanism of the process D + s → a 0 (980) + ρ(770) 0 .However, the interpretation of a 0 (980) + is still controversial [2].Studying the decay D + s → a 0 (980) + ρ(770) 0 can play a key role in understanding the nature of a 0 (980) + .A further motivation to study this decay arises from the tension between theoretical calculations and experimental measurements of the ratio of rates of semileptonic B decays involving leptons from two families.The ratio of the BFs R(D * ) = B(B → D * τ ν)/B(B → D * lν) (l = e, µ) averaged by HFLAV is 0.295 ± 0.011 ± 0.008.It differs from the Standard Model prediction of 0.258 ± 0.005 by 2.6 standard deviations [11].This suggests violation of the lepton flavor universality.However, experimental measurements at LHCb have a large systematic uncertainty from the limited knowledge of the inclusive D + s → π + π + π − X decay [12].A precise measurement of the decay D + s → π + π + π − η may provide useful input to this problem.This Letter reports the first observation of the D + s → π + π + π − η decay, along with an amplitude analysis and a BF measurement of this decay.This analysis uses e + e − collision data samples corresponding to an integrated luminosity of 6.32 fb −1 collected at the center-of-mass energies (E cm ) between 4.178 and 4.226 GeV.Throughout this Letter, charged-conjugate modes and exchange symmetry of two identical π + are always implied.
The BESIII detector and the upgraded multi-gap resistive plate chambers used in the time-of-flight systems are described in Refs.[13,14], respectively.Simulated data samples produced with a geant4-based [15] Monte Carlo (MC) program, which includes the geometric description of the BESIII detector and the detector response, are used to determine detection efficiencies and to estimate backgrounds.The simulation models the beam energy spread and initial state radiation (ISR) in the e + e − annihilations with the generator kkmc [16].The inclusive MC sample includes the production of open charm processes modeled with conexc [17], the ISR production of vector charmonium(-like) states, and the continuum processes incorporated in kkmc [16].The known decay modes are modeled with evtgen [18] using BFs taken from the PDG [2], and the remaining unknown charmonium decays are modeled with lundcharm [19].Final state radiation (FSR) from charged final-state particles is incorporated using photos [20].
We employ the double-tag (DT) technique [21] to select s events.The single-tag (ST) D − s candidates are reconstructed from eight hadronic decay modes (tag side): , and K − π + π − .These tag modes are combined to perform an amplitude analysis.A DT candidate is selected by reconstructing the D + s → π + π + π − η decay (signal side) from the remaining particles that are not used in the ST reconstruction.Here, the K 0 S , π 0 , η and η mesons are reconstructed from K 0 S → π + π − , π 0 → γγ, η → γγ and η → π + π − η decays, respectively.The selection criteria for charged and neutral particle candidates are identical to those used in Ref. [22].For the decay mode D ) is required to be within the range [1.87, 2.06] GeV/c 2 .We calculate the recoiling mass in the e + e − center-of-mass system, where p D − s is the momentum of the reconstructed D − s and m D − s is the known mass of the D − s meson [2].The value of M rec is required to be in the range [2.05, 2.18] GeV/c 2 for the data sample collected at 4.178 GeV to suppress the non-D * ± s D ∓ s events.The M rec ranges for the other data samples are the same as those in Ref. [23].
A seven-constraint kinematic fit is applied to the e + e − → D * ± s D ∓ s → γD + s D − s candidates, where D − s de-cays to one of the tag modes and D + s decays to the signal mode.In addition to the constraints of four-momentum conservation in the e + e − center-of-mass system, the invariant masses of η, tag D − s , and D * + s candidates are constrained to their individual PDG values [2].If there are multiple candidates (in ≈ 15% of the selected events) in an event, the one with the minimum χ 2 of the sevenconstraint kinematic fit is accepted.
To suppress the background from D + s → K 0 S π + η decays, a secondary vertex fit [24] is performed on the π + π − pair.If its invariant mass is in the range [0.487, 0.511] GeV/c 2 and the flight distance between the interaction point (IP) [24] and the decay point is two times greater than its uncertainty, we reject these candidates.To suppress the background from the decay To suppress the background where individual photons from random π 0 's feed into the η → γγ reconstruction, we define two invariant masses M (γ η γ π 0 ) and M (γ η γ other ), where γ η , γ π 0 , and γ other denote the photon of η from the signal side, the photon of π 0 from the tag side, and the other photons including the transition photon from D * +(−) s , respectively.We reject events with M (γ η γ π 0 ) or M (γ η γ other ) in the range [0.115, 0.150] GeV/c 2 .
We further reduce the background by using a gradientboosted decision tree (BDT) implemented in the TMVA software package [25].The BDT takes four discriminating variables as inputs: the recoiling mass of D * s , the momentum of the lower-energy photon from η, the invariant mass of the two photons used to reconstruct η, and the energy of the transition photon from D * s .We place a requirement on the output of the BDT to ensure the samples have a purity greater than 85%: (87.0 ± 3.8)%, (85.6 ± 4.9)%, and (89.2 ± 9.2)% at E cm = 4.178, 4.189-4.219,and 4.226 GeV, respectively.
An unbinned maximum likelihood method is adopted in the amplitude analysis of the D + s → π + π + π − η decay.The likelihood function is constructed with a probability density function (PDF) in which the momenta of the four final-state particles are used as inputs.The total likelihood is the product of the likelihoods for all data samples.The total amplitude is modeled as a coherent sum over all intermediate processes M (p j ) = ρ n e iφn A n (p j ), where ρ n e iφn is the coefficient of the n th amplitude with magnitude ρ n and phase φ n .The n th amplitude A n (p j ) is given by , where the indices 1, 2 and 3 correspond to the two subsequent intermediate resonances and the D + s meson, F i n is the Blatt-Weisskopf barrier factor [26,27] and P i n is the propagator of the intermediate resonance.The function S n describes angular momentum conservation in the decay and is constructed using the covariant tensor formalism [27].The relativistic Breit-Wigner (RBW) [28] function is used to describe the propagator for the resonances η(1405), f 1 (1420) and a 1 (1260).The resonance ρ(770) + is parameterized by the Gounaris-Sakurai [29] lineshape, and the resonances a 0 (980) and f 0 (980) are parameterized by a coupled Flatté formula, and the parameters are fixed to the values given in Ref. [30] and Ref. [31], respectively.We use the same parameterization to describe f 0 (500) as Ref. [32].The masses and widths of the intermediate resonances, except for a 0 (980), f 0 (980) and f 0 (500), are taken from Ref. [2].
The background PDF B(p j ) is constructed from inclusive MC samples by using a multi-dimensional kernel density estimator [33] with five independent Lorentz invariant variables (M π + π + , M π + π − , M π + η , M π − η and M π + π − η ).As a consequence, the combined PDF can be written as where is the acceptance function determined with phase-space (PHSP) MC samples generated with a uniform distribution of the D + s → π + π + π − η decay over PHSP, B is B/ , and R 4 dp j is the element of four-body PHSP.The normalization integral in the denominator is determined by a MC technique as described in Ref. [34].
In the initial amplitude fit, we include a few obvious components.Then, further amplitudes are added one at a time to the fit, and the statistical significance of the new amplitude is calculated with the change of the log-likelihood, after taking the change of the degrees of freedom into account.Only amplitudes with significance larger than 5 σ are chosen for the optimal set.The dominant CF amplitude to this final state is D + s → a 1 (1260) + η, a 1 (1260) → [ρ(770) 0 π + ] S , where the subscript S means that the angular momentum of ρ 0 π + combination is zero (S-wave).Thus, we choose this amplitude as the reference and its phase is fixed to 0. The D + s → a 0 (980) + ρ(770) 0 , a 0 (980) + → π + η decay is observed with a significance larger than 7σ.We also consider some possible amplitudes involving the resonances f 0 (500), f 0 (980), f 1 (1285), η(1295), η(1405), η(1475), f 1 (1420), f 1 (1510) and π(1300), as well as nonresonant components.Moreover, charge conjugation for D + s → η(1405)(f 1 (1420))π + with η(1405)(f 1 (1420)) → a 0 (980) + π − and η(1405)(f 1 (1420)) → a 0 (980) − π + requires their magnitudes and phases to be the same.Finally, eleven amplitudes are retained in the nominal fit, as listed in Table I.The mass projections of the fit are shown in Fig. 2. For the n th component, its contribution relative to the total BF is quantified by the fit fraction (FF) defined by The final amplitudes, their phases and FFs are listed in Table I.The sum of the FFs of all the components is (95.0 ± 4.9)%.A χ 2 value is calculated to quantify the quality of the fit, as defined in Ref. [34].The goodness of fit is χ 2 /NDOF = 153.2/133= 1.15 and the p-value is 11.1%.In addition, 300 sets of signal MC samples with the same size as the data samples are generated to validate the fit performance, as described in Ref. [22].No significant bias is found in our fit.
projections for data with the best fit superimposed.The dots with error bars, the red histograms and the blue histograms are data from all samples, the total fit and the background contribution estimated with inclusive MC samples, respectively.Since the two π + 's in the final state are identical particles, they have a symmetric wave function in some states.The projections of M π + π − , M π + π − η and M π + η are added for the two π + 's.
We determine the systematic uncertainties by taking the differences between the values of φ n and FF n found by the nominal fit and those found from fit variations.The masses and widths of the intermediate states are varied by ±1 σ [2].The masses and coupling constants of a 0 (980) and f 0 (980) are varied within the uncertainties reported in Ref. [30] and Ref. [31], respectively.The barrier radii for D + s and the other intermediate states are varied by ±1 GeV −1 .The uncertainties related to lineshape are estimated by using an alternative RBW function for f 0 (500) with mass and width fixed to 526 MeV/c 2 and 535 MeV [35].The uncertainties from detector effects are investigated with the same method as described in Ref. [22].The uncertainty related to the background is estimated by varying the background yield within statistical uncertainty, and constructing the background PDF with the other five independent variables.The total uncertainties are obtained by adding all the contributions in quadrature, and are listed in Table I.
The systematic uncertainties for the BF measurement are described next.The uncertainty of the signal yield and total ST yield is assigned to be 1.4% by examining the changes of the fit yields when varying the signal and background shapes.The π ± tracking (PID) efficiencies are studied using samples of e + e − → K + K − π + π − (e + e − → K + K − π + π − and π + π − π + π − ) events.The corresponding systematic uncertainties are estimated as 0.3% (0.4%).The uncertainty due to the η reconstruction efficiency is 2.0% [7].The uncertainty from the amplitude model is estimated to be 0.4%, which is the change of signal efficiency when the parameters are varied according to the covariance matrix in the nominal amplitude fit.The uncertainty due to MC simulation sample size is 0.3%, and that from the BF of B(η → γγ) is 0.5% [2].Adding these uncertainties in quadrature gives a total systematic uncertainty of 2.9%.

FIG. 2 :
FIG. 2: The (a)M π + π + , (b) M π + π − , (c) M π + η , (d) M π − η , (e) M π + π + π − and (f) M π + π − ηprojections for data with the best fit superimposed.The dots with error bars, the red histograms and the blue histograms are data from all samples, the total fit and the background contribution estimated with inclusive MC samples, respectively.Since the two π + 's in the final state are identical particles, they have a symmetric wave function in some states.The projections of M π + π − , M π + π − η and M π + η are added for the two π + 's.

FIG. 3 :
FIG. 3: Fits to (a)-(h) the Mtag distributions of the ST candidates and (i) the Msig distribution of the signal candidates.The dots with error bars are data from all samples.The red solid lines are the fits.The blue dashed lines are the fitted background shapes.The pair of pink arrows mark the chosen signal region.In the plots for D − s → K 0 S K − and D − s → π − η decays, the blue dashed lines include contributions from D − → K 0 S π − and D − s → π + π − π − η backgrounds, respectively.

TABLE I :
Phases and FFs for various intermediate processes.The first and the second uncertainties are statistical and systematic, respectively.
These fits result in a total ST yield of Y tag = 479, 093± 1952 and a signal yield of Y sig = 2139 ± 78 events.An updated inclusive MC sample based on our amplitude analysis results is used to determine the ST efficiencies ( i ST ) and DT efficiencies ( i DT ).Inserting these numbers in B(D + s