Observation of $D^+\to\eta\eta\pi^+$ and improved measurement of $D^{0(+)}\to\eta\pi^+\pi^{-(0)}$

Using an $e^+e^-$ annihilation data sample corresponding to an integrated luminosity of $2.93\,\rm fb^{-1}$ collected at the center-of-mass energy of 3.773\,GeV with the BESIII detector, we measure the absolute branching fractions of $D^+\to\eta\eta\pi^+$, $D^+\to\eta\pi^+\pi^0$, and $D^0\to\eta\pi^+\pi^-$ to be $(2.96 \pm 0.24 \pm 0.13)\times 10^{-3}$, $(2.23 \pm 0.15 \pm 0.11)\times 10^{-3}$, and $(1.20 \pm 0.07 \pm 0.04)\times 10^{-3}$, respectively, where the first uncertainties are statistical and the second ones systematic. The $D^+\to\eta\eta\pi^+$ decay is observed for the first time and the branching fractions of $D^{+(0)}\to\eta\pi^+\pi^{0(-)}$ are measured with much improved precision. In addition we test for $CP$ asymmetries in the separated charge-conjugate branching fractions; no evidence of $CP$ violation is found.


I. INTRODUCTION
The goal of the experimental studies of hadronic D meson decays is to explore strong and weak interaction effects. Various experiments have measured the branching fractions (BF) of hadronic decays of D mesons [1]. However, measurements of singly Cabibbosuppressed decays to final states containing one or more η mesons are still limited [1]. Recently, the BESIII Collaboration presented measurements of D 0 → ηπ 0 π 0 and D 0 → ηηπ 0 [2]. The isospin-related decay modes D + → ηπ + π 0 and D 0 → ηπ + π − were measured with large uncertainties by the CLEO Collaboration [3] and there is no measurement for D + → ηηπ + . Improved measurements of D + → ηπ + π 0 , D 0 → ηπ + π − and search for D + → ηηπ + will be useful to clarify the gaps between the inclusive and known exclusive D → ηX decay rates. On the other hand, measurements of these decays provide important inputs for charm and B physics. For instance, these multibody hadronic D decays are crucial backgrounds in semi-tauonic decays of B mesons, thus precision measurements of these hadronic decays are important for the test of lepton flavor universality [4].
This paper presents the first measurement of the branching fraction of D + → ηηπ + and the improved measurements of D +(0) → ηπ + π 0(−) using an e + e − data sample of 2.93 fb −1 taken at the center-of-mass energy √ s = 3.773 GeV [5]. In order to search for CP violation in D decays [6,7], the asymmetries of the BFs of the charge-conjugate decays, defined as A CP = B(D→f )−B(D→f ) B(D→f )+B(D→f ) , have also been measured for the first time. Throughout the paper, charge conjugate modes are implied, except for the A CP measurements.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [8] located at the Beijing Electron Positron Collider (BEPCII) [9]. The cylindrical core of the BESIII detector consists of a helium-based multi-layer 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 muon chambers interleaved with steel. The acceptance of charged particles and photons is 93% of the 4π solid angle. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for the 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 of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps.
Simulated samples produced with the geant4based [10] Monte Carlo (MC) package which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds.
The MC sample used includes production of DD pairs with consideration of quantum coherence for all neutral D modes, the non-DD decays of the ψ(3770), the ISR production of the J/ψ and ψ(3686) states, and the continuum processes incorporated in kkmc [11]. The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modelled with the generator kkmc [11].
The known decay modes of the D mesons and the charmonium states are modelled with evtgen [12] using BFs taken from the Particle Data Group [1], and the remaining unknown decays from the charmonium states with lundcharm [13]. Final state radiation is incorporated with the photos package [14].

III. MEASUREMENT METHOD
Using e + e − annihilations at √ s = 3.773 GeV, we produce DD pairs with no additional hadrons. Events where oneD meson is fully reconstructed are referred to as "single tagged" (ST). A correct tag guarantees the presence of the other D meson, and we search for the hadronic decays D 0(+) → ηπ + π −(0) and D + → ηηπ + recoiling against a taggedD meson. Events with both a tag and such a signal-mode candidate are referred to as "double-tag" events (DT). In this analysis, the tagged D 0 mesons are reconstructed using three hadronic decay modes: K + π − , K + π − π 0 , and K + π − π − π + , while the tagged D − mesons are reconstructed using six hadronic decay modes: K + π − π − , K 0 S π − , K + π − π − π 0 , K 0 S π − π 0 , K 0 S π + π − π − , and K + K − π − . For a specific tag mode i, the yields of the taggedD mesons (N i ST ) and of the DT events (N i DT ) are where N DD is the number of DD pairs, B i ST and B sig are the BFs of theD tag decay mode i and the D signal decay mode, ǫ i ST is the efficiency for finding the tag candidate, ǫ i DT is the efficiency for simultaneously finding the taḡ D and the signal decay. Finally, B sub is the appropriate BF product of η → γγ and π 0 → γγ in the signal decay; i.e., B sub is equal to B 2 η→γγ , B η→γγ B π 0 →γγ , and B η→γγ for D + → ηηπ + , D + → ηπ + π 0 , and D 0 → ηπ + π − , respectively. Combining the above equation, the BF for the signal decay is given by where N ST and N DT are the total ST and DT yields, and ǫ sig is the average efficiency of reconstructing the signal decay (with a tag present), weighted by the measured yields of tag modes in data which is given by (3)

IV. EVENT SELECTION
The event selection criteria used in this work are the same as those used in Refs. [15][16][17][18]. All charged tracks are required to be within a polar-angle (θ) range of |cosθ| < 0.93. Except for those from K 0 S decays, all tracks must originate from an interaction region defined by V xy < 1 cm and V z < 10 cm. Here, V xy(z) is the distance of the closest approach of the charged track to the interaction point perpendicular to (along) the beam.
Charged kaons and pions are identified with the information of the TOF and the dE/dx measured by the MDC. Confidence levels for pion and kaon hypotheses (CL π and CL K ) are calculated.
Kaon and pion candidates are required to satisfy CL K > CL π and CL π > CL K , respectively.
The K 0 S mesons are reconstructed in the decay K 0 S → π + π − . Two oppositely charged tracks are required to satisfy V z < 20 cm, but without V xy and PID requirements. The two tracks are constrained to originate from a common vertex, and their invariant mass is required to satisfy |M π + π − −M K 0 S is the nominal mass [1]. The vertex of K 0 S candidate is required to be more than two standard deviations of the vertex resolution away from the interaction point.
The π 0 and η mesons are reconstructed from their decay into two photons. Photon candidates are selected from the list of EMC showers. The shower time is required to be within 700 ns of the event start time. The shower energy is required to be greater than 25 (50) MeV if the crystal with the maximum energy deposit in that cluster is in the barrel (endcap) region [8]. The opening angle between the candidate shower and the nearest charged track must be greater than 10 • . Photon pairs with an invariant mass in the interval 0.115 < M γγ < 0.150 GeV/c 2 (0.515 < M γγ < 0.570 GeV/c 2 ) are accepted as π 0 (η) candidates. To improve resolution, a one-constraint kinematic fit is imposed on each selected photon pair, in which the γγ invariant mass is constrained to the π 0 or η nominal mass [1].
In the selection of the tagged candidates ofD 0 → K + π − , backgrounds from cosmic rays and Bhabha events must be suppressed. First, the two charged tracks must have a TOF time difference less than 5 ns and they must not be consistent with being a muon pair or an electron−positron pair. Second, there must be at least one EMC shower with an energy larger than 50 MeV or at least one additional charged track detected in the MDC [19]. Also, for the D 0 → ηπ + π − candidate events, the invariant mass of the π + π − combination is required to be outside the mass window of |M π + π − − M K 0 S | < 30 MeV/c 2 to reject the backgrounds from the D 0 → K 0 S η decays. The taggedD (signal D) meson is identified by two variables, the energy difference and the beam-constrained mass where the superscript tag (sig) represents the taggedD candidate and signal D candidate, E beam is the beam energy, p tag (sig) and E tag (sig) are the momentum and energy of theD (D) candidate in the rest frame of e + e − system. For each tag (signal) mode, if there are multiple candidates in an event, only the one with the minimum |∆E tag (sig) | is kept. The tag side is required to satisfy ∆E tag ∈ (−55, +40) MeV for the modes containing a π 0 in the final state and ∆E tag ∈ (−25, +25) MeV for the other modes. The signal side is required to satisfy ∆E sig ∈ (−42, +40) MeV, (−68, +52) MeV, and (−40, +38) MeV for D + → ηηπ + , D + → ηπ + π 0 , and D 0 → ηπ + π − , respectively.

V. SINGLE-TAG AND DOUBLE-TAG YIELDS
The ST yields are obtained from maximum likelihood fits to the M tag BC distributions of the accepted tagged D candidates in data, as shown in Fig. 1. In the fits, theD signal is modeled by an MC-simulated shape convolved with a double Gaussian function describing the resolution difference between data and MC simulation. The combinatorial background shape is described by the ARGUS function [20]. The ST yields and the ST efficiencies are summarized in Table I. The total ST yields (N ST ) are 1558195 ± 2113 for D − and 2386575 ± 1928 for D 0 , where the uncertainties are statistical. These yields are slightly different from those reported in Refs. [15][16][17], due to the lack of M BC window requirements.  BKGII, is from events spread along the diagonal, which are mainly from the e + e − → qq processes. The third one, BKGIII, comes from events with both D andD reconstructed incorrectly which spread out the full plot. To extract the DT yield in data, a two-dimensional (2D) unbinned maximum likelihood fit [21] on this distribution is performed. In the fit, the probability density functions  ) in data, where the uncertainties are statistical. The efficiencies do not include the BFs for K 0 S → π + π − and π 0 → γγ.

Tag mode
527193 ± 761 65.60 ± 0.09 K + π − π 0 1138068 ± 1373 37.69 ± 0.04 K + π − π − π + 721314 ± 1120 38.98 ± 0.06 K + π − π − 798935 ± 1011 51.90 ± 0.08 K 0 S π − 93308 ± 329 51.80 ± 0.17 K + π − π − π 0 258044 ± 1036 26.92 ± 0.09 K 0 S π − π 0 221792 ± 1274 28.27 ± 0.10 K 0 S π − π − π + 115532 ± 645 28.60 ± 0.14 K + K − π − 70548 ± 470 42.13 ± 0.25 (PDFs) of the four components mentioned above are constructed as • signal: a(M sig BC , M tag BC ), where g(x; 0, σ) denotes a Gaussian function with mean of zero and standard deviation of σ, c(x; E beam , ξ, 1 2 ) is an ARGUS function defined as Here, A is a normalization constant (independent for the ARGUS functions in the M sig BC and M tag BC directions) and E beam is the endpoint which is fixed at 1.8865 GeV, and G is the fraction of two Gaussians. The There are some peaking backgrounds in M tag BC vs. M sig BC distribution to consider. For the decay D + → ηηπ + , the peaking backgrounds are from a correct tag with an incorrect signal (D + → π + π 0 π 0 ). For the decay D + → ηπ 0 π + , the peaking backgrounds are from a correct tag with an incorrect signal (D + → K 0 L (K 0 S )π + π 0 , K 0 S → π 0 π 0 , or D + → π + π 0 π 0 ). For these peaking backgrounds, the shapes are modeled based on MC simulation and the normalizations are fixed according to the corresponding BFs in PDG [1]. Figure 3 shows the M tag BC and M sig BC projections of the 2D fits to data. From these 2D fits, we obtain the DT yields for individual signal decays (N DT ) in the fitted M tag (sig) BC region (1.8365, 1.8865) GeV/c 2 , as shown in the second column of Table II. For each signal decay mode, the statistical significance is calculated according to −2ln(L 0 /L max ), where L max is the maximal likelihood of the nominal fit and L 0 is the likelihood of the corresponding fit without the signal component. The statistical significance for the three signal decays are all found to be greater than 10σ.

VI. BRANCHING FRACTIONS
To ensure the reliability of signal efficiency, we have examined the M ηP , M ηπ + , and M P π + distributions of D → ηP π + candidate events after requiring |M sig BC − M D | < 0.006 GeV/c 2 . Here, P denotes the daughter particles of η, π 0 , and π − for D + → ηηπ + , D + → ηπ + π 0 , and D 0 → ηπ + π − decays, respectively. Figure 4 shows the Dalitz plots of three signal decay modes in data, and there are no significant ρ 0,± and a 0 (980) 0,± signals in these Dalitz plots. However, due to some possible resonances, the phase-space MC distributions of M ηP , M ηπ + , and M P π + do not agree well with the data distributions. To solve this problem, the MC generator is modified to produce the correct invariant mass distributions according to the Dalitz plots in data. In the Dalitz plot, the background component is modeled by the inclusive MC simulation, while the signal components generated according to an efficiencycorrected MC simulation. These modified MC samples are in good agreement with the data distributions and are therefore used to determine the averaged efficiencies of the signal decays (ǫ sig ), which are summarized in Table II.
The absolute BFs of the signal decays obtained according to Eq. (2), are summarized in Table II Table III.

VII. SYSTEMATIC UNCERTAINTIES
With the DT method, most of uncertainties related to the taggedD are canceled. A summary of the systematic uncertainties in the BF measurements is given in Table IV and are discussed below.   [3]. The efficiencies do not include the BFs of η → γγ and π 0 → γγ. The uncertainties in NDT and ǫsig are statistical. The first and second uncertainties of Bsig and BCLEO are statistical and systematic, respectively. 24.96 ± 0.12 2.96 ± 0.24 ± 0.13 N/A D + → ηπ + π 0 381 ± 26 28.11 ± 0.13 2.23 ± 0.15 ± 0.11 1.38 ± 0.31 ± 0.16 D 0 → ηπ + π − 450 ± 25 39.98 ± 0.17 1.20 ± 0.07 ± 0.04 1.09 ± 0.13 ± 0.09 TABLE III: Summary of the ST yields (NST), the signal yields (NDT), and the signal efficiencies (ǫsig) used to determinate the BFs (Bsig) and CP asymmetries (ACP ) for D → sig andD → sig. For ACP , the first and second uncertainties are statistical and systematic, respectively. The uncertainties for other values are only statistical.  • ST yields: The uncertainties in the total ST yields come from the fits to the M BC spectra of the tagged D 0 and D − candidates. They have been previously estimated to be 0.5% for both neutral and charged D in Refs. [15][16][17].
• Tracking (PID) of π ± : The tracking (PID) efficiencies of π ± are investigated with DT DD hadronic events by using a partial reconstruction technique. The systematic uncertainty for each charged particle due to tracking (PID) is estimated to be 0.5%.
The momentum weighted data-MC difference in π 0 reconstruction efficiencies is found to be (−0.5 ± 1.0)%, where the uncertainty is statistical. After correcting the MC efficiencies by the momentum weighted data-MC difference in π 0 reconstruction efficiency, the systematic uncertainty due to π 0 reconstruction is assigned as 1.0% per π 0 . The systematic uncertainty due to η reconstruction is assumed to be the same as π 0 reconstruction and fully correlated.
• 2D yield fits: The systematic uncertainty due to the 2D fits of the M tag BC vs. M sig BC distributions is evaluated by repeating the measurements with an alternative fit range of (1.8300, 1.8865) GeV/c 2 , an alternative signal shape with different MC matching requirements, alternative endpoints of the ARGUS function, E beam ± 0.2 MeV/c 2 , and with the quoted BFs of peaking backgrounds varied by ±1σ. The total systematic uncertainties are assigned based on the changes of the BFs from each of these sources summed in quadrature, yielding 1.0%, 2.1%, and 0.8% for D + → ηηπ + , D + → ηπ + π 0 , and D 0 → ηπ + π − , respectively.
• ∆E sig requirement: The systematic uncertainties due to the ∆E sig requirement are assigned by comparing the DT efficiencies with and without smearing by the data-MC difference of the ∆E sig resolution for the signal MC events. Here, the ∆E sig resolution differences are obtained by using larger DT samples of D 0 → K − π + η, D 0 → K 0 S η, and D + → π + π 0 π 0 with the same tags. The maximum change of the DT efficiency is taken to be the systematic uncertainties, which is 0.3% for D + → ηηπ + , D + → ηπ + π 0 , and D 0 → ηπ + π − .
• Modified MC generator: The systematic uncertainty in the modified MC generator is studied with an alternative input Dalitz plot obtained by varying the MC-simulated background sizes.
• K 0 S rejection: The efficiency uncertainty from K 0 S rejection is estimated by using an alternative rejection window of ±40 MeV/c 2 around the K 0 S nominal mass. The change in the BF, 1.4%, is assigned as the systematic uncertainty for D 0 → ηπ + π − .
• Asymmetry of CP ± components: The measurement of the BF of D 0 → ηπ + π − is affected by CP ± eigenstate components in the D 0 → ηπ + π − decay. The asymmetry of CP + and CP − components in this decay is examined by the CP + tag of D 0 → K + K − and the CP − tag of D 0 → K 0 S π 0 . Combined with the strong-phase factors of the flavor tagsD 0 → K − π + ,D 0 → K − π + π 0 , and D 0 → K − π + π + π − [1,22,23], the impact on the BF of D 0 → ηπ + π − is found to be (1.0±0.9)% with the same method described in Ref. [24]. After correcting the BF of D 0 → ηπ + π − by this factor, 0.9% is assigned as an associated uncertainty.
• Measurement method: The reliability of the measurement method has been validated with ten sets of inclusive MC samples. Each data set has equivalent integrated luminosity of data. It is found that the measured BFs of D + → ηηπ + , D + → ηπ + π 0 , and D 0 → ηπ + π − are shifted by 3.0%, 1.8%, and 1.7%, respectively, which are directly taken as a systematic uncertainty due to measurement method.
The total systematic uncertainty obtained by adding the above contributions in quadrature is 4.5%, 4.9%, and 3.7% for D + → ηηπ + , D + → ηπ + π 0 , and D 0 → ηπ + π − , respectively. In the determinations of A CP , the uncertainties of π 0 and η reconstruction, quoted BFs, MC modeling, measurement method for each decay, π + π − tracking and PID as well as strong phase for D 0 (D 0 ) → ηπ + π − are assumed to cancel, while for D +/− → ηπ +/− π 0 and ηηπ +/− decays, the uncertainties of π +/− tracking and PID are assumed to be un-canceled. The remaining systematic uncertainties have been estimated separately with the same methods mentioned above. With current statistics, no evidence of CP violation is found.

VIII. SUMMARY
With a data sample corresponding to an integrated luminosity of 2.93 fb −1 taken at √ s = 3.773 GeV with the BESIII detector, we measure the absolute BFs of the singly Cabibbo-suppressed decays D + → ηηπ + , D + → ηπ + π 0 , and D 0 → ηπ + π − . The BF of D + → ηηπ + is measured for the first time. The BFs of D + → ηπ + π 0 and D 0 → ηπ + π − are consistent with the CLEOc's results [3] within 2.2σ and 0.6σ, respectively. The asymmetries of the BFs of D andD decays in the three channels have also been examined, no evidence of CP violation is found. In the near future, amplitude analyses of these three decays with larger data samples at BESIII and Belle II will offer opportunity to explore two-body decays D → ρη, a 0 (980)π, and a 0 (980)η.