Measurements of the absolute branching fractions of $D^{0(+)}\to K\bar{K}\pi\pi$ decays

Multi-body hadronic $D^{0(+)}$ decays provide an ideal laboratory to study strong and weak interactions. Amplitude analyses of these decays offer comprehensive information of quasi-two-body $D^{0(+)}$ decays, which are important to explore $D\bar{D}^{0}$ mixing, charge-parity ($CP$) violation and quark SU(3)-flavor asymmetry breaking phenomenon ref5 ; theory_1 ; theory_2 ; chenghy1 ; yufs . In particular, for the search of $CP$ violation, it is important to understand the intermediate structures for the singly Cabibbo-suppressed decays of $D^{0(+)}\to K\bar{K}\pi\pi$ xwkang ; Charles:2009ig ; yufs-cpv .

Current measurements of the $D^{0(+)}\to K\bar{K}\pi\pi$ decays containing $K^{0}_{S}$ or $\pi^{0}$ are limited pdg2018 . The branching fractions (BFs) of $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{-}$ FOCUS_kskspipi ; ARGUS_kkpipi , $D^{+}\to K^{0}_{S}K^{-}\pi^{+}\pi^{+}$ FOCUS_kskpipi , $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$ FOCUS_kskpipi , and $D^{+}\to K^{+}K^{-}\pi^{+}\pi^{0}$ ACCMOR_kkpipi0 were only determined relative to some well known decays or via topological normalization, with poor precision. This paper presents the first direct measurements of the absolute BFs for the decays $D^{0}\to K^{+}K^{-}\pi^{0}\pi^{0}$, $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{-}$, $D^{0}\to K^{0}_{S}K^{-}\pi^{+}\pi^{0}$, $D^{0}\to K^{0}_{S}K^{+}\pi^{-}\pi^{0}$, $D^{+}\to K^{+}K^{-}\pi^{+}\pi^{0}$, $D^{+}\to K^{0}_{S}K^{+}\pi^{0}\pi^{0}$, $D^{+}\to K^{0}_{S}K^{-}\pi^{+}\pi^{+}$, $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$, and $D^{+}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{0}$. The $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{0}\pi^{0}$ decay is not included since it suffers from poor statistics and high background. Throughout this paper, charge conjugate processes are implied. An $e^{+}e^{-}$ collision data sample corresponding to an integrated luminosity of 2.93 fb^-1 lum_bes3 collected at a center-of-mass energy of $\sqrt{s}=$ 3.773 GeV with the BESIII detector is used to perform this analysis.

The BESIII detector is a magnetic spectrometer BESIII located at the Beijing Electron Positron Collider (BEPCII) Yu:IPAC2016-TUYA01 . The cylindrical core of the BESIII detector 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 identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over $4\pi$ solid angle. The charged-particle momentum resolution at $1~{}{\rm 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 geant4-based geant4 Monte Carlo (MC) package including the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds. The simulation includes the beam-energy spread and initial-state radiation (ISR) in the $e^{+}e^{-}$ annihilations modeled with the generator kkmc kkmc . The inclusive MC samples consist of the production of $D\bar{D}$ pairs with consideration of quantum coherence for all neutral $D$ modes, the non-$D\bar{D}$ decays of the $\psi(3770)$, the ISR production of the $J/\psi$ and $\psi(3686)$ states, and the continuum processes. The known decay modes are modeled with evtgen evtgen using the BFs taken from the Particle Data Group (PDG) pdg2018 , and the remaining unknown decays from the charmonium states are modeled with lundcharm lundcharm . The final-state radiations from charged final-state particles are incorporated with the photos package photos .

The $D^{0}\bar{D}^{0}$ or $D^{+}D^{-}$ pair is produced without an additional hadron in $e^{+}e^{-}$ annihilations at $\sqrt{s}=3.773$ GeV. This process offers a clean environment to measure the BFs of the hadronic $D$ decay with the double-tag (DT) method. The single-tag (ST) candidate events are selected by reconstructing a $\bar{D}^{0}$ or $D^{-}$ in the following hadronic final states: $\bar{D}^{0}\to K^{+}\pi^{-}$, $K^{+}\pi^{-}\pi^{0}$, and $K^{+}\pi^{-}\pi^{-}\pi^{+}$, and $D^{-}\to K^{+}\pi^{-}\pi^{-}$, $K^{0}_{S}\pi^{-}$, $K^{+}\pi^{-}\pi^{-}\pi^{0}$, $K^{0}_{S}\pi^{-}\pi^{0}$, $K^{0}_{S}\pi^{+}\pi^{-}\pi^{-}$, and $K^{+}K^{-}\pi^{-}$. The event in which a signal candidate is selected in the presence of an ST $\bar{D}$ meson, is called a DT event. The BF of the signal decay is determined by

where $N^{\rm tot}_{\rm ST}=\sum_{i}N_{{\rm ST}}^{i}$ and $N^{\rm net}_{\rm DT}$ are the total yields of the ST and DT candidates in data, respectively. $N_{{\rm ST}}^{i}$ is the ST yield for the tag mode $i$. For the signal decays involving $K^{0}_{S}$ meson(s) in the final states, $N^{\rm net}_{\rm DT}$ is the net DT yields after removing the peaking background from the corresponding non-$K^{0}_{S}$ decays. For the other signal decays, the variable corresponds to the fitted DT yields as described later. Here, $\epsilon_{{\rm sig}}$ is the efficiency of detecting the signal $D$ decay, averaged over the tag mode $i$, which is given by:

where $\epsilon^{i}_{{\rm ST}}$ and $\epsilon^{i}_{{\rm DT}}$ are the efficiencies of detecting ST and DT candidates in the tag mode $i$, respectively.

The selection criteria of $K^{\pm}$, $\pi^{\pm}$, $K^{0}_{S}$, and $\pi^{0}$ are the same as those used in the analyses presented in Refs. epjc76 ; cpc40 ; bes3-pimuv ; bes3-Dp-K1ev ; bes3-etaetapi ; bes3-omegamuv ; bes3-etamuv ; bes3-etaX . All charged tracks, except those from $K^{0}_{S}$ decays, are required to have a polar angle $\theta$ with respect to the beam direction within the MDC acceptance $|\rm{cos\theta}|<0.93$, and a distance of closest approach to the interaction point (IP) within 10 cm along the beam direction and within 1 cm in the plane transverse to the beam direction. Particle identification (PID) for charged pions, kaons, and protons is performed by exploiting TOF information and the specific ionization energy loss $dE/dx$ measured by the MDC. The confidence levels for pion and kaon hypotheses ($CL_{\pi}$ and $CL_{K}$) are calculated. Kaon and pion candidates are required to satisfy $CL_{K}>CL_{\pi}$ and $CL_{\pi}>CL_{K}$, respectively.

The $K^{0}_{S}$ candidates are reconstructed from two oppositely charged tracks to which no PID criteria are applied and which masses are assumed to be that of pions. The charged tracks from the $K^{0}_{S}$ candidate must satisfy $|\rm{cos\theta}|<0.93$. In addition, due to the long lifetime of the $K^{0}_{S}$ meson, there is a less stringent criterion on the distance of closest approach to the IP in the beam direction of less than 20 cm and no requirement on the distance of closest approach in the plane transverse to the beam direction. Furthermore, the $\pi^{+}\pi^{-}$ pairs are constrained to originate from a common vertex and their invariant mass is required to be within $(0.486,0.510)~{}{\rm GeV}/c^{2}$, which corresponds to about three times the fitted resolution around the nominal $K^{0}_{S}$ mass. The decay length of the $K^{0}_{S}$ candidate is required to be larger than two standard deviations of the vertex resolution away from the IP.

The $\pi^{0}$ candidate is reconstructed via its $\gamma\gamma$ decay. The photon candidates are selected using the information from the EMC shower. It is required that each EMC shower starts within 700 ns of the event start time and its energy is greater than 25 (50) MeV in the barrel (end cap) region of the EMC BESIII . The energy deposited in the nearby TOF counters is included to improve the reconstruction efficiency and energy resolution. The opening angle between the candidate shower and the nearest charged track must be greater than $10^{\circ}$. The $\gamma\gamma$ pair is taken as a $\pi^{0}$ candidate if its invariant mass is within $(0.115,\,0.150)$ GeV$/c^{2}$. To improve the resolution, a kinematic fit constraining the $\gamma\gamma$ invariant mass to the $\pi^{0}$ nominal mass pdg2018 is imposed on the selected photon pair.

To select $\bar{D}^{0}\to K^{+}\pi^{-}$ candidates, the backgrounds from cosmic rays and Bhabha events are rejected by using the same requirements described in Ref. deltakpi . In the selection of $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{-}\pi^{+}$ candidates, the $\bar{D}^{0}\to K^{0}_{S}K^{\pm}\pi^{\mp}$ decays are suppressed by requiring the mass of all $\pi^{+}\pi^{-}$ pairs to be outside $(0.478,0.518)$ GeV/$c^{2}$.

The tagged $\bar{D}$ mesons are identified using two variables, namely the energy difference

and the beam-constrained mass

Here, $E_{\rm b}$ is the beam energy, $\vec{p}_{\rm tag}$ and $E_{\rm tag}$ are the momentum and energy of the $\bar{D}$ candidate in the rest frame of $e^{+}e^{-}$ system, respectively. For each tag mode, if there are multiple candidates in an event, only the one with the smallest $|\Delta E_{\rm tag}|$ is kept. The tagged $\bar{D}$ candidates are required to satisfy $\Delta E_{\rm tag}\in(-55,40)$ MeV for the tag modes containing $\pi^{0}$ in the final states and $\Delta E_{\rm tag}\in(-25,25)$ MeV for the other tag modes, thereby taking into account the different resolutions.

To extract the yields of ST $\bar{D}$ mesons for individual tag modes, binned-maximum likelihood fits are performed on the corresponding $M_{\rm BC}^{\rm tag}$ distributions of the accepted ST candidates following Refs. epjc76 ; cpc40 ; bes3-pimuv ; bes3-Dp-K1ev ; bes3-etaetapi ; bes3-omegamuv ; bes3-etamuv . In the fits, the $\bar{D}$ 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 an ARGUS function ARGUS defined as $c_{f}(f;E_{\rm end},\xi_{f})=A_{f}\cdot f\cdot\sqrt{1-\frac{f^{2}}{E^{2}_{\rm end}}}\cdot\exp\left[\xi_{f}\left(1-\frac{f^{2}}{E^{2}_{\rm end}}\right)\right]$, where $f$ denotes $M^{\rm tag}_{\rm BC}$, $E_{\rm end}$ is an endpoint fixed at 1.8865 GeV, $A_{f}$ is a normalization factor, and $\xi_{f}$ is a free parameter. The resulting fits to the $M_{\rm BC}$ distributions for each mode are shown in Fig. 1. The total yields of the ST $\bar{D}^{0}$ and $D^{-}$ mesons in data are $2327839\pm 1860$ and $1558159\pm 2113$, respectively, where the uncertainties are statistical only.

In the recoiling sides against the tagged $\bar{D}$ candidates, the signal $D$ decays are selected by using the residual tracks that have not been used to reconstruct the tagged $\bar{D}$ candidates. To suppress the $K^{0}_{S}$ contribution in the individual mass spectra for the $D^{0}\to K^{+}K^{-}\pi^{0}\pi^{0}$, $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{-}$, and $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$ decays, the $\pi^{+}\pi^{-}$ and $\pi^{0}\pi^{0}$ invariant masses are required to be outside $(0.468,0.528)$ GeV/$c^{2}$ and $(0.438,0.538)$ GeV/$c^{2}$, respectively. To suppress the background from $D^{0}\to K^{-}\pi^{+}\omega$ in the identification of the $D^{0}\to K^{0}_{S}K^{-}\pi^{+}\pi^{0}$ process, the $K^{0}_{S}\pi^{0}$ invariant mass is required to be outside $(0.742,0.822)$ GeV/$c^{2}$. These requirements correspond to at least five times the fitted mass resolution away from the fitted mean of the mass peak.

The signal $D$ mesons are identified using the energy difference $\Delta E_{\rm sig}$ and the beam-constrained mass $M_{\rm BC}^{\rm sig}$, which are calculated with Eqs. (3) and (4) by substituting “tag” with “sig”. For each signal mode, if there are multiple candidates in an event, only the one with the smallest $|\Delta E_{\rm sig}|$ is kept. The signal decays are required to satisfy the mode-dependent $\Delta E_{\rm sig}$ requirements, as shown in the second column of Table 1. To suppress incorrectly identified $D\bar{D}$ candidates, the opening angle between the tagged $\bar{D}$ and the signal $D$ is required to be greater than $160^{\circ}$, resulting in a loss of (2-6)% of the signal and suppressing (8-55)% of the background.

Figure 2 shows the $M_{\rm BC}^{\rm tag}$ versus $M_{\rm BC}^{\rm sig}$ distribution of the accepted DT candidates in data. The signal events concentrate around $M_{\rm BC}^{\rm tag}=M_{\rm BC}^{\rm sig}=M_{D}$, where $M_{D}$ is the nominal $D$ mass pdg2018 . The events with correctly reconstructed $D$ ($\bar{D}$) and incorrectly reconstructed $\bar{D}$ ($D$), named BKGI, are spread along the lines around $M_{\rm BC}^{\rm tag}=M_{D}$ or $M_{\rm BC}^{\rm sig}=M_{D}$. The events smeared along the diagonal, named BKGII, are mainly from the $e^{+}e^{-}\to q\bar{q}$ processes. The events with uncorrelated and incorrectly reconstructed $D$ and $\bar{D}$, named BKGIII, disperse across the whole allowed kinematic region.

For each signal $D$ decay mode, the yield of DT events ($N^{\rm fit}_{\rm DT}$) is obtained from a two-dimensional (2D) unbinned maximum-likelihood fit cleo-2Dfit on the $M_{\rm BC}^{\rm tag}$ versus $M_{\rm BC}^{\rm sig}$ distribution of the accepted candidates. In the fit, the probability density functions (PDFs) of signal, BKGI, BKGII, and BKGIII are constructed as

• •

  signal: $a(x,y)$,

• •

  BKGI: $b(x)\cdot c_{y}(y;E_{\rm b},\xi_{y})+b(y)\cdot c_{x}(x;E_{\rm b},\xi_{x})$,

• •

  BKGII: $c_{z}(z;\sqrt{2}E_{\rm b},\xi_{z})\cdot g(k)$, and

• •

  BKGIII: $c_{x}(x;E_{\rm b},\xi_{x})\cdot c_{y}(y;E_{\rm b},\xi_{y})$,

respectively. Here, $x=M_{\rm BC}^{\rm sig}$, $y=M_{\rm BC}^{\rm tag}$, $z=(x+y)/\sqrt{2}$, and $k=(x-y)/\sqrt{2}$. The PDFs of signal $a(x,y)$, $b(x)$, and $b(y)$ are described by the corresponding MC-simulated shapes. $c_{f}(f;E_{\rm end},\xi_{f})$ is an ARGUS function ARGUS defined above, where $f$ denotes $x$, $y$, or $z$; $E_{\rm b}$ is fixed at 1.8865 GeV. $g(k)$ is a Gaussian function with mean of zero and standard deviation parametrized by $\sigma_{k}=\sigma_{0}\cdot(\sqrt{2}E_{\rm b}/c^{2}-z)^{p}$, where $\sigma_{0}$ and $p$ are fit parameters.

Combinatorial $\pi^{+}\pi^{-}$ pairs from the decays $D^{0}\to K^{0}_{S}2(\pi^{+}\pi^{-})$ [and $D^{0}\to 3(\pi^{+}\pi^{-})$], $D^{0}\to K^{-}\pi^{+}\pi^{+}\pi^{-}\pi^{0}$, $D^{0}\to K^{+}\pi^{+}\pi^{-}\pi^{-}\pi^{0}$, $D^{+}\to K^{-}\pi^{+}\pi^{+}\pi^{+}\pi^{-}$, $D^{+}\to K^{+}2(\pi^{+}\pi^{-})$, $D^{+}\to K^{+}\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, $D^{+}\to K^{0}_{S}\pi^{+}\pi^{+}\pi^{-}\pi^{0}$ [and $D^{+}\to 2(\pi^{+}\pi^{-})\pi^{+}\pi^{0}$] may also satisfy the $K^{0}_{S}$ selection criteria and form peaking backgrounds around $M_{D}$ in the $M_{\rm BC}^{\rm sig}$ distributions for $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{-}$, $D^{0}\to K^{0}_{S}K^{-}\pi^{+}\pi^{0}$, $D^{0}\to K^{0}_{S}K^{+}\pi^{-}\pi^{0}$, $D^{+}\to K^{0}_{S}K^{+}\pi^{0}\pi^{0}$ $D^{+}\to K^{0}_{S}K^{-}\pi^{+}\pi^{+}$, $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$, and $D^{+}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{0}$, respectively. This kind of peaking background is estimated by selecting events in the $K^{0}_{S}$ sideband region of $(0.454,0.478)\cup(0.518,0.542)~{}{\rm GeV}/c^{2}$. For $D^{0}\to K^{0}_{S}K^{-}\pi^{+}\pi^{0}$, $D^{0}\to K^{0}_{S}K^{+}\pi^{-}\pi^{0}$, $D^{+}\to K^{0}_{S}K^{-}\pi^{+}\pi^{+}$, $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$, and $D^{+}\to K^{0}_{S}K^{+}\pi^{0}\pi^{0}$ decays, one-dimensional (1D) signal and sideband regions are used. For $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{-}$ and $D^{+}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{0}$ decays, 2D signal and sideband regions are used. The 2D $K^{0}_{S}$ signal region is defined as the square region with both $\pi^{+}\pi^{-}$ combinations lying in the $K^{0}_{S}$ signal regions. The 2D $K^{0}_{S}$ sideband 1 (2) regions are defined as the square regions with 1 (2) $\pi^{+}\pi^{-}$ combination(s) located in the 1D $K^{0}_{S}$ sideband regions and the rest in the 1D $K^{0}_{S}$ signal region. Figure 3 shows 1D and 2D $\pi^{+}\pi^{-}$ invariant-mass distributions as well as the $K^{0}_{S}$ signal and sideband regions.

For the signal decays involving $K^{0}_{S}$ meson(s) in the final states, the net yields of DT events are calculated by subtracting the sideband contribution from the DT fitted yield by

Here, $N=1$ for the decays with one $K^{0}_{S}$ meson while $N=2$ for the decays with two $K^{0}_{S}$ mesons. The combinatorial $\pi^{+}\pi^{-}$ backgrounds are assumed to be uniformly distributed and double-counting is avoided by subtracting (2) yields from (1) yields appropriately. $N^{\rm fit}_{\rm DT}$ and $N^{\rm fit}_{{\rm sid}i}$ are the fitted $D$ yields in the 1D or 2D signal region and sideband $i$ region, respectively. For the other signal decays, the net yields of DT events are $N^{\rm fit}_{\rm DT}$. Figure 4 shows the $M^{\rm tag}_{\rm BC}$ and $M^{\rm sig}_{\rm BC}$ projections of the 2D fits to data. From these 2D fits, we obtain the DT yields for individual signal decays as shown in Table 1.

For each signal decay mode, the statistical significance is calculated according to $\sqrt{-2{\rm ln({\mathcal{L}_{0}}/{\mathcal{L}_{\rm max}}})}$, where ${\mathcal{L}}_{\rm max}$ and ${\mathcal{L}}_{0}$ are the maximum likelihoods of the fits with and without involving the signal component, respectively. The effect of combinatorial $\pi^{+}\pi^{-}$ backgrounds in the $K^{0}_{S}$-signal regions has been considered for the decays involving a $K^{0}_{S}$. The statistical significance for each signal decay is found to be greater than $8\sigma$.

Each of the $D^{0}\to K^{0}_{S}K^{-}\pi^{+}\pi^{0}$, $D^{+}\to K^{+}K^{-}\pi^{+}\pi^{0}$, $D^{+}\to K^{0}_{S}K^{-}\pi^{+}\pi^{+}$, and $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$ decays is modeled by the corresponding mixed signal MC samples, in which the dominant decay modes containing resonances of $K^{*}(892)$, $\rho(770)$, and $\phi$ are mixed with the phase space (PHSP) signal MC samples. The mixing ratios are determined by examining the corresponding invariant mass and momentum spectra. The other decays, which are limited in statistics, are generated with the PHSP generator. The momentum and the polar angle distributions of the daughter particles and the invariant masses of each two- and three-body particle combinations of the data agree with those of the MC simulations. As an example, Fig. 5 shows the invariant mass distributions of two or three-body particle combinations of $D^{+}\to K^{+}K^{-}\pi^{+}\pi^{0}$ candidate events for data and MC simulations.

The measured values of $N^{\rm net}_{{\rm DT}}$, $\epsilon_{{\rm sig}}$, and the obtained BFs are summarized in Table 1. The current world-average values are also given for comparison. The signal efficiencies have been corrected by the necessary data-MC differences in the selection efficiencies of $K^{\pm}$ and $\pi^{\pm}$ tracking and PID procedures and the $\pi^{0}$ reconstruction. These efficiencies also include the BFs of the $K^{0}_{S}$ and $\pi^{0}$ decays. The efficiency for $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$ ($D^{0}\to K^{0}_{S}K^{-}\pi^{+}\pi^{0}$) is lower than that of $D^{+}\to K^{0}_{S}K^{-}\pi^{+}\pi^{+}$ ($D^{0}\to K^{0}_{S}K^{+}\pi^{-}\pi^{0}$) due to the $K^{0}_{S}$ $(\omega)$ rejection in the $\pi^{+}\pi^{-}$ ($K^{0}_{S}\pi^{0}$) mass spectrum.

The systematic uncertainties are estimated relative to the measured BFs and are discussed below. In BF determinations using Eq. (1), all uncertainties associated with the selection of tagged $\bar{D}$ canceled in the ratio. The systematic uncertainties in the total yields of ST $\bar{D}$ mesons related to the $M_{\rm BC}$ fits to the ST $\bar{D}$ candidates, were previously estimated to be 0.5% for both neutral and charged $\bar{D}$ epjc76 ; cpc40 ; bes3-pimuv .

The tracking and PID efficiencies for $K^{\pm}$ or $\pi^{\pm}$, $\epsilon_{K\,{\rm or}\,\pi}^{\rm tracking\,(PID)}[{\rm data}]$ and $\epsilon_{K\,{\rm or}\,\pi}^{\rm tracking\,(PID)}[{\rm MC}]$, are investigated using DT $D\bar{D}$ hadronic events. The averaged ratios between data and MC efficiencies ($f_{K\,{\rm or}\,\pi}^{\rm tracking\,(PID)}=\epsilon_{K\,{\rm or}\,\pi}^{\rm tracking\,(PID)}[{\rm data}]/\epsilon_{K\,{\rm or}\,\pi}^{\rm tracking\,(PID)}[{\rm MC}]$) of tracking (PID) for $K^{\pm}$ or $\pi^{\pm}$ are weighted by the corresponding momentum spectra of signal MC events, giving $f_{K}^{\rm tracking}$ to be $1.022{\text{-}}1.031$ and $f_{\pi}^{\rm tracking}$ to be close to unity. After correcting the MC efficiencies by $f_{K}^{\rm tracking}$, the residual uncertainties of $f_{K\,{\rm or}\,\pi}^{\rm tracking}$ are assigned as the systematic uncertainties of tracking efficiencies, which are (0.4-0.7)% per $K^{\pm}$ and (0.2-0.3)% per $\pi^{\pm}$. $f_{K}^{\rm PID}$ and $f_{\pi}^{\rm PID}$ are all close to unity and their individual uncertainties, (0.2-0.3)%, are taken as the associated systematic uncertainties per $K^{\pm}$ or $\pi^{\pm}$.

The systematic error related to the uncertainty in the $K_{S}^{0}$ reconstruction efficiency is estimated from measurements of $J/\psi\to K^{*}(892)^{\mp}K^{\pm}$ and $J/\psi\to\phi K_{S}^{0}K^{\pm}\pi^{\mp}$ control samples sysks and found to be 1.6% per $K^{0}_{S}$. The systematic uncertainty of $\pi^{0}$ reconstruction efficiency is assigned as (0.7-0.8)% per $\pi^{0}$ from a study of DT $D\bar{D}$ hadronic decays of $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{0}$ and $\bar{D}^{0}\to K^{0}_{S}\pi^{0}$ decays tagged by either $D^{0}\to K^{-}\pi^{+}$ or $D^{0}\to K^{-}\pi^{+}\pi^{+}\pi^{-}$ epjc76 ; cpc40 .

The systematic uncertainty in the 2D fit to the $M_{\rm BC}^{\rm tag}$ versus $M_{\rm BC}^{\rm sig}$ distribution is examined via the repeated measurements in which the signal shape and the endpoint of the ARGUS function ($\pm 0.2$ MeV/$c^{2}$) are varied. Quadratically summing the changes of the BFs for these two sources yields the corresponding systematic uncertainties.

The systematic uncertainty due to the $\Delta E_{\rm sig}$ requirement is assigned to be 0.3%, which corresponds to the largest efficiency difference with and without smearing the data-MC Gaussian resolution of $\Delta E_{\rm sig}$ for signal MC events. Here, the smeared Gaussian parameters are obtained by using the samples of DT events $D^{0}\to K^{0}_{S}\pi^{0}$, $D^{0}\to K^{-}\pi^{+}\pi^{0}$, $D^{0}\to K^{-}\pi^{+}\pi^{0}\pi^{0}$, and $D^{+}\to K^{-}\pi^{+}\pi^{+}\pi^{0}$ versus the same $\bar{D}$ tags in our nominal analysis. The systematic uncertainties due to $K^{0}_{S}$ sideband choice and $K^{0}_{S}$ rejection mass window are assigned by examining the changes of the BFs via varying nominal $K^{0}_{S}$ sideband and corresponding rejection window by $\pm 5$ MeV/$c^{2}$.

For the decays whose efficiencies are estimated with mixed signal MC events, the systematic uncertainty in the MC modeling is determined by comparing the signal efficiency when changing the percentage of MC sample components. For the decays whose efficiencies are estimated with PHSP-distributed signal MC events, the uncertainties are assigned as the change of the signal efficiency after adding the possible decays containing $K^{*}(892)$ or $\rho(770)$. The imperfect simulations of the momentum and $\cos\theta$ distributions of charged particles are considered as a source of systematic uncertainty. The signal efficiencies are re-weighted by those distributions in data with background subtracted. The largest change of the re-weighted to nominal efficiencies, 0.9%, is assigned as the corresponding systematic uncertainty.

The measurements of the BFs of the neutral $D$ decays are affected by quantum correlation effect. For each neutral $D$ decay, the $CP$-even component is estimated by the $CP$-even tag $D^{0}\to K^{+}K^{-}$ and the $CP$-odd tag $D^{0}\to K^{0}_{S}\pi^{0}$. Using the same method as described in Ref. QC-factor and the necessary parameters quoted from Refs. R-ref1 ; R-ref2 ; R-ref3 , we find the correction factors to account for the quantum correlation effect on the measured BFs are $(98.3^{+1.6}_{-1.1{\,\rm stat}})\%$, $(98.1^{+2.8}_{-1.7{\,\rm stat}})\%$, $(95.9^{+3.4}_{-2.7{\,\rm stat}})\%$, and $(98.4^{+1.1}_{-1.0{\,\rm stat}})\%$ for $D^{0}\to K^{+}K^{-}\pi^{0}\pi^{0}$, $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{-}$, $D^{0}\to K^{0}_{S}K^{-}\pi^{+}\pi^{0}$, and $D^{0}\to K^{0}_{S}K^{+}\pi^{-}\pi^{0}$, respectively. After correcting the signal efficiencies by the individual factors, the residual uncertainties are assigned as systematic uncertainties.

The uncertainties due to the limited MC statistics for various signal decays, (0.4-0.8)%, are taken into account as a systematic uncertainty. The uncertainties of the quoted BFs of the $K^{0}_{S}\to\pi^{+}\pi^{-}$ and $\pi^{0}\to\gamma\gamma$ decays are 0.07% and 0.03%, respectively pdg2018 .

The efficiencies of $D\bar{D}$ opening angle requirement is studied by using the DT events of $D^{0}\to K^{-}\pi^{+}\pi^{+}\pi^{-}$, $D^{0}\to K^{-}\pi^{+}\pi^{0}\pi^{0}$, and $D^{+}\to K^{-}\pi^{+}\pi^{+}\pi^{0}$ tagged by the same tag modes in our nominal analysis. The difference of the accepted efficiencies between data and MC simulations, 0.4% for the decays without $\pi^{0}$, 0.8% for the decays involving one $\pi^{0}$ and 0.3% for the decays involving two $\pi^{0}$s, is assigned as the associated systematic uncertainty.

Table 2 summarizes the systematic uncertainties in the BF measurements. For each signal decay, the total systematic uncertainty is obtained by adding the above effects in quadrature to be (2.6-6.0)% for various signal decay modes.

In summary, by analyzing a data sample obtained in $e^{+}e^{-}$ collisions at $\sqrt{s}=3.773$ GeV with the BESIII detector and corresponding to an integrated luminosity of 2.93 fb^-1, we obtained the first direct measurements of the absolute BFs of nine $D^{0(+)}\to K\bar{K}\pi\pi$ decays containing $K^{0}_{S}$ or $\pi^{0}$ mesons. The $D^{0}\to K^{+}K^{-}\pi^{0}\pi^{0}$, $D^{0}\to K^{0}_{S}K^{-}\pi^{+}\pi^{0}$, $D^{0}\to K^{0}_{S}K^{+}\pi^{-}\pi^{0}$, $D^{+}\to K^{0}_{S}K^{+}\pi^{0}\pi^{0}$, and $D^{+}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{0}$ decays are observed for the first time. Compared to the world-average values, the BFs of the $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{-}$, $D^{+}\to K^{+}K^{-}\pi^{+}\pi^{0}$, $D^{+}\to K^{0}_{S}K^{-}\pi^{+}\pi^{+}$, and $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$ decays are measured with improved precision. Our BFs of $D^{+}\to K^{0}_{S}K^{-}\pi^{+}\pi^{+}$ and $D^{+}\to K^{0}_{S}K^{+}\pi^{+}\pi^{-}$ are in agreement with individual world averages within $1\sigma$ while our BFs of $D^{0}\to K^{0}_{S}K^{0}_{S}\pi^{+}\pi^{-}$ and $D^{+}\to K^{+}K^{-}\pi^{+}\pi^{0}$ deviate with individual world averages by $2.3\sigma$ and $2.8\sigma$, respectively. The precision of the BF of $D^{+}\to K^{+}K^{-}\pi^{+}\pi^{0}$ is improved by a factor of about seven. Future amplitude analyses of all these $D^{0(+)}\to K\bar{K}\pi\pi$ decays with larger data samples foreseen at BESIII bes3-white-paper , Belle II belle2-white-paper , and LHCb lhcb-white-paper will supply rich information of the two-body decay modes containing scalar, vector, axial and tensor mesons, thereby benefiting the understanding of quark SU(3)-flavor symmetry.

Authors thank for valuable discussions with Prof. Fu-sheng Yu. The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support. This work is supported in part by National Key Basic Research Program of China under Contract No. 2015CB856700; National Natural Science Foundation of China (NSFC) under Contracts Nos. 11775230, 11475123, 11625523, 11635010, 11735014, 11822506, 11835012, 11935015, 11935016, 11935018, 11961141012; the Chinese Academy of Sciences (CAS) Large-Scale Scientific Facility Program; Joint Large-Scale Scientific Facility Funds of the NSFC and CAS under Contracts Nos. U1532101, U1932102, U1732263, U1832207; CAS Key Research Program of Frontier Sciences under Contracts Nos. QYZDJ-SSW-SLH003, QYZDJ-SSW-SLH040; 100 Talents Program of CAS; INPAC and Shanghai Key Laboratory for Particle Physics and Cosmology; ERC under Contract No. 758462; German Research Foundation DFG under Contracts Nos. Collaborative Research Center CRC 1044, FOR 2359; Istituto Nazionale di Fisica Nucleare, Italy; Ministry of Development of Turkey under Contract No. DPT2006K-120470; National Science and Technology fund; STFC (United Kingdom); The Knut and Alice Wallenberg Foundation (Sweden) under Contract No. 2016.0157; The Royal Society, UK under Contracts Nos. DH140054, DH160214; The Swedish Research Council; U. S. Department of Energy under Contracts Nos. DE-FG02-05ER41374, DE-SC-0012069.