Measurement of branching fractions and search for $CP$ violation in $D^{0}\to\pi^{+}\pi^{-}\eta$, $D^{0}\to K^{+}K^{-}\eta$, and $D^{0}\to\phi\eta$ at Belle

We measure the branching fractions and $CP$ asymmetries for the singly Cabibbo-suppressed decays $D^{0}\to\pi^{+}\pi^{-}\eta$, $D^{0}\to K^{+}K^{-}\eta$, and $D^{0}\to\phi\eta$, using 980 fb$^{-1}$ of data from the Belle experiment at the KEKB $e^+e^-$ collider. We obtain \begin{eqnarray} \mathcal{B}(D^{0}\to\pi^{+}\pi^{-}\eta)&=&[1.22\pm 0.02\,({\rm stat})\pm 0.02\,({\rm syst})\pm 0.03\,(\mathcal{B}_{\rm ref})]\times 10^{-3}\,, \nonumber \\ \mathcal{B}(D^{0}\to K^{+}K^{-}\eta)&=&[1.80\,^{+0.07}_{-0.06}\,({\rm stat})\pm 0.04\,({\rm syst})\pm 0.05\,(\mathcal{B}_{\rm ref})]\times 10^{-4}\,, \nonumber \\ \mathcal{B}(D^{0}\to\phi\eta)&=&[1.84\pm 0.09\,({\rm stat})\pm 0.06\,({\rm syst})\pm 0.05\,(\mathcal{B}_{\rm ref})]\times 10^{-4}\,, \nonumber \end{eqnarray} where the third uncertainty ($\mathcal{B}_{\rm ref}$) is from the uncertainty in the branching fraction of the reference mode $D^{0}\to K^{-}\pi^{+}\eta$. The color-suppressed decay $D^{0}\to\phi\eta$ is observed for the first time, with very high significance. The results for the $CP$ asymmetries are \begin{eqnarray} A_{CP}(D^{0}\to\pi^{+}\pi^{-}\eta)&=&[0.9\pm 1.2\,({\rm stat})\pm 0.5\,({\rm syst})]\%\,, \nonumber \\ A_{CP}(D^{0}\to K^{+}K^{-}\eta)&=&[-1.4\pm 3.3\,({\rm stat})\pm 1.1\,({\rm syst})]\%\,, \nonumber \\ A_{CP}(D^{0}\to\phi\eta)&=&[-1.9\pm 4.4\,({\rm stat})\pm 0.6\,({\rm syst})]\%\,. \nonumber \end{eqnarray} The results for $D^{0}\to\pi^{+}\pi^{-}\eta$ are a significant improvement over previous results. The branching fraction and $A_{CP}$ results for $D^{0}\to K^{+}K^{-}\eta$, and the $A_{CP}$ result for $D^{0}\to\phi\eta$, are the first such measurements. No evidence for $CP$ violation is found in any of these decays.


Introduction
Singly Cabibbo-suppressed (SCS) decays of charmed mesons provide a promising opportunity to study CP violation in the charm sector. Within the Standard Model, CP violation in charm decays is expected to be of the order of 10 −3 or smaller [1,2], and thus challenging to observe. SCS decays are of special interest, as interference that includes a new physics amplitude could lead to large CP violation. The CP asymmetry between D 0 → f and D 0 →f decays (A CP ) is defined as The only observation of CP violation in the charm sector to date is from the LHCb experiment, where a difference in A CP between the SCS D 0 → K + K − and D 0 → π + π − decays [3] was observed: ∆A CP = (−15.4 ± 2.9) × 10 −4 . In this paper, we investigate two analogous SCS decays, D 0 → π + π − η and D 0 → K + K − η. A search for a CP asymmetry in the first decay was performed by the BESIII experiment; the resulting precision was 6% [4]. There have been no results for D 0 → K + K − η decays to date. Theoretically, it is difficult to predict CP asymmetries for three-body decays, while some predictions exist for intermediate two-body processes: A CP (D 0 → ρ 0 η) is predicted to be −0.53 × 10 −3 from tree amplitudes alone, and −0.23 × 10 −3 after considering QCD-penguin and weak penguin-annihilation [1]. The asymmetry A CP (D 0 → φη) is predicted to be zero in several theoretical models [1]. A precise measurement of branching fractions (B) for these three-body decays is an important step towards searching for CP violation in these channels.
To identify the flavor of the neutral D meson when produced, we reconstruct D * + → D 0 π + s and D * − → D 0 π − s decays; the charge of the daughter π ± s (which has low momentum and is referred to as the "slow" pion) identifies whether the D meson is D 0 or D 0 . The raw asymmetry measured (A raw ) receives contributions from several sources: where A D 0 →f CP is the CP asymmetry for D 0 → f ; A D * + FB is the forward-backward asymmetry due to γ-Z 0 interference and higher-order QED effects [10] in e + e − → cc collisions; and A πs ε is the asymmetry resulting from a difference in reconstruction efficiencies between π + s and π − s . This asymmetry depends on the transverse momentum p T (π s ) and polar angle θ(π s ) of the π s in the laboratory frame. We correct for this by weighting signal events by a factor [1 − A πs ε (p T , cos θ)] for D 0 decays, and by a factor [1 + A πs ε (p T , cos θ)] for D 0 decays. After this weighting, we are left with the π s -corrected asymmetry A corr (cos θ * ) = A CP + A FB (cos θ * ) . (

1.3)
Since A FB is an odd function of the cosine of the D * + polar angle θ * in the e + e − centerof-mass (CM) frame, and A CP is independent of cos θ * , we extract A CP and A FB (cos θ * ) via

4)
A FB (cos θ * ) = A corr (cos θ * ) − A corr (− cos θ * ) 2 . (1.5) Fitting the values of A CP for different cos θ * bins to a constant gives our final measurement of A CP for D 0 → f .

Belle detector and data sets
This measurement is based on the full data set of the Belle experiment, which corresponds to a total integrated luminosity of 980 fb −1 [11] collected at or near the Υ(nS) (n = 1, 2, 3, 4, 5) resonances. The Belle experiment ran at the KEKB energy-asymmetric collider [12,13]. The Belle detector is a large-solid-angle magnetic spectrometer consisting 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 comprising CsI(Tl) crystals located inside a superconducting solenoid coil providing a 1.5 T magnetic field. An iron flux-return located outside the coil is instrumented to detect K 0 L mesons and to identify muons. A detailed description of the detector is given in Refs. [11,14].
We use Monte Carlo (MC) simulated events to optimize selection criteria, study backgrounds, and evaluate the signal reconstruction efficiency. Signal MC events are generated by EVTGEN [15] and propagated through a detector simulation based on GEANT3 [16]. Final-state radiation from charged particles is simulated using the PHOTOS package [17]. Three-body decays are generated according to phase space. An MC sample of "generic" events, corresponding to an integrated luminosity four times that of the data, is used to develop selection criteria. It includes BB events and continuum processes e + e − → qq, where q = u, d, s, c. At the Υ(5S) resonance, the MC includes B is the kinetic energy released in the D * + decay.

Event selection and optimization
We reconstruct the signal decays D 0 → π + π − η and D 0 → K + K − η, and the reference decay D 0 → K − π + η, in which the D 0 originates from D * + → D 0 π + , as follows. 1 Charged tracks are identified as K ± or π ± candidates using a likelihood ratio R K ≡ L K /(L K + L π ), where L K (L π ) is the likelihood that a track is a K ± (π ± ) based on the photon yield in the ACC, dE/dx information in the CDC, and time-of-flight information from the TOF [18]. Tracks having R K > 0.60 are identified as K ± candidates; otherwise, they are considered as π ± candidates. The corresponding efficiencies are approximately 90% for kaons and 95% for pions. Tracks that are highly electron-like (R e > 0.95) or muon-like (R µ > 0.95) are rejected, where the electron and muon likelihood ratios R e and R µ are determined mainly using information from the ECL and KLM detectors, respectively [19,20]. Charged tracks are required to have at least two SVD hits in the +z direction (defined as the direction opposite that of the positron beam), and at least two SVD hits in the x-y (transverse) plane. The nearest approach of the π + s track to the e + e − interaction point (IP) is required to be less than 1.0 cm in the x-y plane, and less than 3.0 cm along the z axis.
Photon candidates are identified as energy clusters in the ECL that are not associated with any charged track. The photon energy (E γ ) is required to be greater than 50 MeV in the barrel region (covering the polar angle 32 • < θ < 129 • ), and greater than 100 MeV in the endcap region (12 • < θ < 31 • or 132 • < θ < 157 • ). The ratio of the energy deposited in the 3 × 3 array of crystals centered on the crystal with the highest energy, to the energy deposited in the corresponding 5 × 5 array of crystals, is required to be greater than 0.80.
Candidate η → γγ decays are reconstructed from photon pairs having an invariant mass satisfying 500 MeV/c 2 < M (γγ) < 580 MeV/c 2 . This range corresponds to about 3σ in M (γγ) resolution. The absolute value of the cosine of the η → γ 1 γ 2 decay angle, defined as cos , is required to be less than 0.85. This retains around 89% of the signal while reducing backgrounds by a factor of two. To further suppress backgrounds, we remove η candidates in which both photon daughters can be combined with other photons in the event to form π 0 → γγ candidate decays satisfying |M γγ − m π 0 | < 10 MeV/c 2 , where m π 0 is the nominal π 0 mass [5]. This veto requirement has an efficiency of 95% while reducing backgrounds by a factor of three (D 0 → K + K − η) and four (D 0 → π + π − η).
Candidate D 0 → π + π − η, D 0 → K + K − η, and D 0 → K − π + η decays are reconstructed by combining π ± and K ± tracks with η candidates. A vertex fit is performed with the two charged tracks to obtain the D 0 decay vertex position; the resulting fit quality is labeled χ 2 v . To improve the momentum resolution of the η, the γ daughters are subjected to a fit in which the photons are required to originate from the D 0 vertex position, and the invariant mass is constrained to be that of the η meson [5]. The fit quality of this mass constraint (χ 2 m ) is required to satisfy χ 2 m < 8, and the resulting η momentum is required to be greater than 0.70 GeV/c. For D 0 → π + π − η candidates, we veto events [5], to suppress background from CF D 0 → K 0 S η decays. This veto range corresponds to about 3σ in resolution. The D 0 invariant mass M is required to satisfy 1.850 GeV/c 2 < M < 1.878 GeV/c 2 for D 0 → K + K − η candidates; 1.840 GeV/c 2 < M < 1.884 GeV/c 2 for D 0 → π + π − η candidates; and 1.842 GeV/c 2 < M < 1.882 GeV/c 2 for D 0 → K − π + η candidates. These ranges correspond to about 2σ in resolution.
Candidate D * + → D 0 π + s decays are reconstructed by combining D 0 candidates with π + s tracks. We first fit for D * + decay vertex using the D 0 momentum vector and decay vertex position, and the IP as a constraint (i.e., the D * + nominally originates from the IP). The resulting goodness-of-fit is labeled χ 2 IP . To improve the resolution in Q, another vertex fit is performed: in this case we constrain the π + s daughter to originate from the D * + decay vertex, and the resulting fit quality is labeled χ 2 s . The sum of the above three fit qualities, , is required to be less than 50; this requirement has a signal efficiency of about 97%. Those D * + candidates satisfying 0 < Q < 15 MeV are retained for further analysis. To eliminate D * + candidates originating from B decays, and to also suppress combinatorial background, the D * + momentum in the CM frame is required to be greater than 2.70 GeV/c.
After the above selection criteria are applied, about 2.1% of D 0 → π + π − η events, 1.3% of D 0 → K − π + η events, and < 0.1% of D 0 → K + K − η events have two or more D * + candidates. For such multi-candidate events, we choose a single candidate: that which has the smallest value of the sum χ 2 vtx + χ 2 m (η). This criterion, according to MC simulation, identifies the correct candidate 54% of the time.

Measurement of the branching fractions
We extract the signal yield via an unbinned maximum-likelihood fit to the Q distribution. The probability density function (PDF) used for signal events is taken to be the sum of a bifurcated Student's t-function (S bif ), which is defined in appendix A, and one or two asymmetric Gaussians (G asym ), with all having a common mean. The PDF used for D 0 → K + K − η signal events is simply an S bif function. These PDFs are explicitly In these expressions, δ i is an asymmetry parameter characterizing the difference between left-side and right-side widths: σ R,L = σ(1±δ). Most of these parameters are fixed to values obtained from MC simulation. However, the parameters µ, σ 0 , and, for the higher-statistics D 0 → K − π + η channel, n l,h , are floated to account for possible differences in resolution between data and MC. For backgrounds, the PDF is taken to be a threshold function f (Q) = Q α e −βQ ; for the CF mode D 0 → K − π + η, we include an additional symmetric Gaussian to describe a small background component originating from misreconstructed D 0 → K − π + π 0 π 0 decays. The parameters of this Gaussian are fixed to values obtained from MC simulation, while all other parameters are floated. No other peaking backgrounds, such as misreconstructed D 0 decays or signal decays in which a pion from the D 0 is swapped with that from the D * + decay, are found in the MC simulation.
The results of the fit are shown in figure 1, along with the pull (N data − N fit )/σ, where σ is the error on N data . All fit residuals look satisfactory. The signal yields in the fitted region 0 < Q < 15 MeV, and in the signal region |Q − 5.86| < 0.80 MeV, are listed in table 1.

Region
Component To measure the branching fraction, we must divide these signal yields by their reconstruction efficiencies. However, the reconstruction efficiency for a decay can vary across the Dalitz plot of three-body phase space, and the Dalitz-plot distribution of D 0 → π + π − η and D 0 → K + K − η decays has not been previously measured. Thus, to avoid systematic uncertainty due to the unknown Dalitz distribution (or decay model), we correct our signal -5 -  yields for reconstruction efficiencies as follows. We divide the Dalitz plot of the data into bins of M 2 (h + h − ) and M 2 (h − η), where h = π or K, determine the reconstruction efficiency independently for each bin, and calculate the corrected signal yield via the formula where i runs over all bins. The values of The formula (4.4) has the following terms: • ε i is the signal reconstruction efficiency for bin i, as determined from a large sample of MC events. The resolutions in Q of the MC samples are adjusted to match those of the data. The efficiencies for D 0 → π + π − η are plotted in figure 2(a), and those for D 0 → K + K − η are plotted in figure 3(a). These efficiencies include a small (∼ 2%) correction for K ± and π ± particle identification (PID) efficiencies, to account for small differences observed between data and MC simulation. This correction is determined using a sample of D * + → [D 0 → K − π + ]π + s decays.
• N tot i is the number of events in the Q signal region and the i th bin of the Dalitz plot. These yields are plotted in figure 2(b) for D 0 → π + π − η and in figure 3 • N bkg is the total background yield in the Q signal region, as obtained from fitting the Q distribution (see figure 1).
• f bkg i is the fraction of background in the i th -bin, with i f i = 1. These fractions are obtained from the Dalitz plot distribution of events in the Q sideband region 2.5 MeV < |Q − 5.86| < 4.9 MeV. The distribution of sideband events is shown in figure 2 There are 10 × 10 = 100 bins in total for D 0 → π + π − η, and 5 × 5 = 25 bins total for D 0 → K + K − η. The final corrected yields obtained using eq. (4.4)  π + π − , M 2 π − η , and M 2 π + η , respectively. Points with error bars show events in the signal region; blue-filled histograms show the estimated background (see text).
The branching fraction of a signal mode relative to that of the normalization mode is determined from the ratio of their respective efficiency-corrected yields: where h = K or π. The efficiency-corrected yield for the normalization channel D 0 → K − π + η is evaluated in a different manner than those of the signal modes. As the Dalitz plot of D 0 → K − π + η decays has been measured with high statistics [6], we use the resulting decay model to generate an MC sample, and use that sample to evaluate the overall reconstruction efficiency. The result, including the small PID efficiency correction, is ε Kπη = (6.870 ± 0.014)%. Dividing the fitted yield for D 0 → K − π + η (see table 1) by this value gives N cor (D 0 → K − π + η) = (2.365 ± 0.011) × 10 6 .
Inserting all efficiency-corrected yields into eq. (4.5) gives the ratios of branching fractions The second error listed is the systematic uncertainty, which is evaluated below (section 4.3).
For D 0 → K + K − η events, the M 2 K + K − distribution shows the D 0 → φη, φ → K + K − decay process to be dominant. However, a non-φ contribution is also visible. We thus The signal yield is obtained as before by fitting the Q distribution. The result is 599 ± 45 events in the signal region, as shown in figure 1(d). The change in likelihood, with and without including a signal component in such Q fitting, is ∆ ln L = 214. As the number of degrees of freedom for the fit with no signal component is three less than that for the nominal fit (parameters N sig , µ, and σ 0 are dropped), this value of ∆ ln L corresponds to a statistical significance for the signal of 20σ.
We divide the signal yields obtained for bins of the Dalitz plot by the efficiencies for these bins [see eq. (4.4)] to obtain N cor = 12443 +1071 −893 (for D 0 → K + K − η with the φ excluded). Thus This result is somewhat higher (but more precise) than a similar measurement by BESIII, (0.59 ± 0.19) × 10 −4 [7].

Measurement of B(D 0 → φη)
As shown in figure 3, the decay D 0 → K + K − η is dominated by the Cabibbo-and colorsuppressed decay D 0 → φη, φ → K + K − . We thus measure the branching fraction for D 0 → φη by performing a two-dimensional fit to the M KK and Q distributions of D 0 → K + K − η events. The fitted region is M KK < 1.08 GeV/c 2 and Q < 15 MeV. In this region, signal decays and background are straightforward to identify: the non-φ D 0 → K + K − η component peaks in Q but not in M KK , whereas combinatorial background containing φ → K + K − decays peak in M KK but not in Q. The signal PDF for M KK is taken to be the sum of a Gaussian and two asymmetric Gaussians, with a common mean for the M (φ) peak. The PDF for Q is taken to be a bifurcated Student's t-function: P sig (Q) = S bif (µ q , σ 0 , δ 0 , n l , n h ) , (4.13) (4.14) Most of the signal shape parameters are fixed to MC values; the parameters µ m , σ m0 for M KK and µ q , σ 0 for Q are floated to account for possible differences in resolution between data and MC simulation. The non-φ D 0 → K + K − η component includes several processes such as D 0 → a 0 (980)η, D 0 → K * 0 (1430) ± K ∓ , and non-resonant D 0 → K + K − η decays. We parameterize this component in M KK by a threshold function, and in Q by the same PDF as that used for the signal: where the threshold value m 0 = 2m K = 0.987354 GeV/c 2 . We consider possible interference between the non-φ component and the D 0 → φη signal as a systematic uncertainty. For combinatorial background, the Q distribution is parameterized with a threshold function. The M KK distribution has two parts: (1) a φ-peak, which is taken to be the same as that of signal decay, and (2) a threshold function. These PDFs take the forms  The relative fraction f φ and the Q threshold parameters are floated; all other parameters are fixed to MC values. The results of the two-dimensional likelihood fit are shown in figure 4. We obtain a signal yield N sig = 728 ± 36 in the full fitted region, and N sig = 600 ± 29 in the signal region |Q − 5.86| < 0.8 MeV and |M K + K − − m φ | < 10 MeV/c 2 . The difference in likelihood, with and without including a signal component, is ∆ ln L = 464.8. As the number of degrees of freedom for the fit with no signal component is one less than that for the nominal fit (parameter N sig is dropped), this value of ∆ ln L corresponds to a statistical significance for D 0 → φη of 31σ.
We evaluate the signal reconstruction efficiency using a large MC sample of D 0 → φη decays. We obtain, for events in the M KK -Q signal region, an efficiency ε = (5.262 ± 0.021)%. Thus, N cor (D 0 → φη, φ → K + K − ) = (1.140 ± 0.055) × 10 4 , and   where the systematic uncertainty includes the small uncertainty on B(φ → K + K − ). This result is consistent with, but notably more precise than, the current world average of (1.8 ± 0.5) × 10 −4 [5]. It is also consistent with theoretical predictions [22,23]. As a consistency check, we calculate the branching fraction of the non-φ D 0 → K + K − η component by subtracting the D 0 → φη branching fraction from the total D 0 → K + K − η result: which is very close to our measurement of B(D 0 → K + K − η) φ−excluded in eq. (4.11).

Systematic uncertainties
The sources of systematic uncertainty in measuring the branching fractions are listed in table 2. These uncertainties are evaluated as follows.
Systematic sources • A correction for PID efficiency is applied to K ± and π ± tracks, to account for a difference in efficiency between data and MC simulation. The correction depends on track momentum and is small; the uncertainty on the correction is even smaller, in the range (0.90-0.97)%. When evaluating this uncertainty for a ratio of branching fractions, we conservatively assume the efficiency corrections for K + and π + tracks (which appear separately in numerator and denominator, or vice-versa) are anticorrelated.
• The uncertainty due to the parameters fixed in the fit for the signal yield is evaluated as follows. We sample these parameters simultaneously from Gaussian distributions, accounting for their correlations, and re-fit for the signal yield. The procedure is repeated 1000 times and these yields are plotted. The ratio of the root-mean-square (RMS) to the mean value of the resulting distribution of signal yields is taken as the systematic uncertainty due to the fixed parameters. The Gaussian distributions from which the parameters are sampled have mean values equal to the fixed values of the parameters, and widths equal to their respective uncertainties.
• For background PDFs, all parameters are floated in the fits except for those describing the amount of background and its shape for D 0 → K − π + η, which are taken from MC simulation. We evaluate this uncertainty due to these fixed parameters as done above, by simultaneously sampling these parameters from Gaussian distributions having mean values equal to the fixed values and widths equal to their respective uncertainties. The RMS of the resulting distribution of D 0 → K − π + η yields is taken as the systematic uncertainty due to the peaking background.
• We correct for differences in mass resolutions (including M , Q and M KK ) between data and MC when calculating reconstruction efficiencies (in eq. (4.4) for D 0 → π + π − η and D 0 → K + K − η, as well as for D 0 → φη and the reference mode). We take the systematic uncertainty of this procedure to be the difference in the ratio of efficiency-corrected signal yields to that of the reference mode obtained both with and without this resolution correction.
• The efficiency for D 0 → π + π − η and D 0 → K + K − η decays is evaluated in bins of the Dalitz plot; see eq. (4.4). This efficiency has uncertainty arising from the number of bins used, from the efficiency values ε i for the various bins, and from the bin-by-bin background subtraction.
to evaluate the effect of uncertainties in ε i , we sample the ε i from Gaussian distributions having mean values equal to the nominal values, and widths equal to their uncertainties. For each sampling, we re-calculate the yield N cor using eq. (4.4) and plot the result. The RMS of this distribution is taken as the systematic uncertainty: 0.21% for D 0 → π + π − η, and 0.43% for D 0 → K + K − η.
the background subtraction procedure depends on the distribution of background over the Dalitz plot. We take this distribution from a data Q sideband region. To evaluate the uncertainty due to this Dalitz distribution, we shift the Q sideband region used by ±0.4 MeV and repeat the procedure. The change is assigned as the systematic uncertainty: 0.02% for D 0 → π + π − η, and 0.03% for D 0 → K + K − η.
• The efficiency for D 0 → K − π + η is evaluated from MC simulation using a Dalitz decay model. The uncertainty due to this model is evaluated by varying the model and recalculating ε. Specifically, our nominal model uses eight intermediate resonances as measured in ref. [6]; as an alternative, we include all thirteen intermediate resonances listed in ref. [6]. The resulting change in our reconstruction efficiency is very small (δε/ε < 0.01%). As a cross check, we calculate N cor for D 0 → K − π + η using eq. (4.4); the result is (2.369 ± 0.007) × 10 6 , which is almost identical with our nominal result (the difference is much smaller than the uncertainty).
• There are small uncertainties in the reconstruction efficiencies ε, which are evaluated from MC simulation, due to the finite statistics of the MC samples used.
• There is an uncertainty arising from the K 0 S veto required for D 0 → π + π − η decays. We evaluate this by changing the veto region from m K 0 S ± 10 MeV/c 2 (∼ 3σ) to m K 0 S ± 15 MeV/c 2 (∼ 5σ); the resulting change in the D 0 → π + π − η signal yield is taken as the systematic uncertainty.
• We consider possible interference between the D 0 → φη amplitude and that of non-φ D 0 → K + K − η. Such interference could alter the M KK distribution used to fit for the D 0 → φη yield. We evaluate this effect by introducing a relative phase θ between the two amplitudes, which modifies the PDF to be In this expression, A φ is a relativistic Breit-Wigner function, F non−φ is the shape function used in the nominal fit to model non-φ decays, and k is a factor that adds a phase shift of π depending on whether the cosine of the K + helicity angle (θ h ) is positive or negative, i.e., k = 0 for cos θ h < 0, and k = 1 otherwise [24]. The helicity angle θ h is defined as the angle between the momentum of the K + and the η in the K + K − rest frame. After fitting with this PDF, we calculate the D 0 → φη yield as the product of the total yield obtained and the fraction f φ given by The result is that the D 0 → φη yield decreases by 2.5%, and thus we assign this value as the systematic uncertainty due to possible interference.
The total systematic uncertainty is obtained by adding in quadrature all the above contributions. The results are listed in table 2.
5 Measurement of CP asymmetries To measure the CP asymmetries, we divide the sample for each channel into D 0 and D 0 decays, where the flavor of the D 0 or D 0 is tagged by the charge of the π ± s from the D * + → D 0 π + s or D * − → D 0 π − s decay. To correct for an asymmetry in π ± s reconstruction efficiencies, we weight events according to the π ± s efficiency mapping of ref. [25]. We simultaneously fit the Q distributions of these weighted samples for D 0 and D 0 decays with parameters fixed in the same way as done for the branching fraction measurements. The parameters N sig and A corr are fitted, where the D 0 and D 0 signal yields are given by N sig (D 0 , D 0 ) = (N sig /2) · (1 ± A corr ). The fit results are shown in figure 5. We obtain N sig = 12975 ± 198 and A corr = (1.44 ± 1.24)% for D 0 → π + π − η decays, and N sig = 1482 ± 60 and A corr = (−0.25 ± 2.96)% for D 0 → K + K − η decays.
These values for A corr include the forward-backward asymmetry A FB . We correct for A FB by calculating A corr in eight bins of cos θ * , where θ * is the polar angle of the D * + with respect to the +z axis in the e + e − CM frame. (d) Figure 5. Simultaneous fit for D 0 → π + π − η (a) and D 0 → π + π − η (b) candidates; and D 0 → A FB are plotted in figures 7(a, d) for D 0 → π + π − η and in figures 7(b, e) for D 0 → K + K − η.

Fitting the A CP values to constants yields the final results
A CP (D 0 → π + π − η) = [0.9 ± 1.2 (stat) ± 0.5 (syst)]% , (5.1) The second error listed is the systematic uncertainty, which is evaluated below (section 5.3). The first result is a factor of four more precise than a recent measurement by BESIII [4], while the latter result is the first such measurement. The A FB values plotted in figures 7(de) decrease with cos θ * and are consistent with (somewhat lower in D 0 → π + π − η than) the leading-order prediction [26] at √ s = 10.6 GeV of A cc FB = −0.029 · cos θ * /(1 + cos 2 θ * ), at the current level of statistics.

Measurement of
We repeat the above procedure to determine A CP for D 0 → φη decays. Here, to determine parameters N sig and A corr , we perform a two-dimensional fit in [M KK , Q] for the D 0 and D 0 samples simultaneously. We allow N sig and A corr to float separately for the φη and nonresonant K + K − η components. The projections of the fit result are shown in figure 6, and the results are N sig = 728 ± 36 and A corr = (−0.17 ± 4.44)%. We perform this fit separately to obtain the A corr values for the eight bins of cos θ * and use eqs. ( where the second error listed is the systematic uncertainty, evaluated below (section 5.3). This result is consistent with zero, as expected [22].      (d) Figure 6. Projections in M K + K − (left) and Q (right) of the two-dimensional (M K + K − , Q) fit, for D 0 → φη (top) and D 0 → φη (bottom). In both cases, the M K + K − projection corresponds to the Q signal region (|Q − 5.86| < 0.8 MeV), and the Q projection corresponds to the M K + K − signal region (|M K + K − − m φ | < 0.01 GeV/c 2 ).

Systematic uncertainties
Fortunately, most systematic uncertainties in measuring A CP cancel. The remaining sources of systematic uncertainty are listed in table 3 and are evaluated as follows.
• There is an uncertainty arising from fixed parameters in the fit used to describe signal and background shapes. We evaluate this uncertainty using the sampling method described previously to evaluate uncertainties for the branching fraction measurement. The resulting uncertainties for A CP are small: < 0.001 for both D 0 → π + π − η and Sources D 0 → K + K − η, and 0.002 for D 0 → φη. We also consider different possible Q and M KK resolutions for the D 0 and D 0 samples by allowing the σ 0 parameter in eqs. (4.1 -4.3) and the σ m0 and σ 0 parameters in eqs. (4.12, 4.13) to vary between the two samples. The change in A CP is 0.004 for D 0 → π + π − η, 0.010 for D 0 → K + K − η, and 0.006 for D 0 → φη. Combining these two uncertainties in quadrature gives the values listed in table 3. . The resulting change in A CP is taken as the systematic uncertainty due to this source. There is a small uncertainty arising from a difference in the detector acceptance near the boundaries cos θ * = ±1; we evaluate this by considering only events with | cos θ * | < 0.90.

• We extract
• We correct for a possible asymmetry in π ± s reconstruction efficiencies by weighting events according to a mapping of efficiencies ε(π s ). There are 56 bins in this map, and the efficiencies for each bin has some uncertainty. We thus vary these efficiencies individually by their uncertainties to create 56 new efficiency maps with +1σ shifts and 56 maps with −1σ shifts. We subsequently weight the D 0 and D 0 samples by these efficiency maps and repeat the fit for N sig and A cor CP . The resulting deviations from the nominal fit result are summed in quadrature to give the systematic uncertainty arising from this source. We obtain +0.050 −0.087 % for D 0 → K + K − η, +0.065 −0.072 % for D 0 → K + K − η, and +0.043 −0.100 % for D 0 → φη. The total systematic uncertainty is obtained by adding in quadrature all the above contributions. The results are listed in table 3.

Conclusion
In summary, based on a data set corresponding to an integrated luminosity of 980 fb −1 recorded by the Belle experiment, we report measurements of the branching fractions of the SCS decays D 0 → π + π − η and D 0 → K + K − η relative to that for the CF decay D 0 → K − π + η. We also measure the relative branching fraction for the resonant decay D 0 → φη; this measurement uses an order of magnitude more data than used for our previous measurement [8] and supersedes it. Our results are: where the third uncertainty is due to the branching fraction for the reference mode D 0 → K − π + η. These results are the most precise to date. The time-integrated CP asymmetries are measured to be A CP (D 0 → π + π − η) = [0.9 ± 1.2 (stat) ± 0.4 (syst)]% , (6.7) A CP (D 0 → K + K − η) = [−1.4 ± 3.3 (stat) ± 1.0 (syst)]% , (6.8) A CP (D 0 → φη) = [−1.9 ± 4.4 (stat) ± 0.6 (syst)]% . (6.9) The first result represents a significant improvement in precision over the previous result [4]. The latter two are the first such measurements. No evidence for CP violation is found.