Amplitude Analysis and Branching Fraction Measurement of $D^0\rightarrow K^-\pi^+\pi^0\pi^0$

Utilizing the dataset corresponding to an integrated luminosity of $2.93$ fb$^{-1}$ at $\sqrt{s}=3.773$ GeV collected by the BESIII detector, we report the first amplitude analysis and branching fraction measurement of the $D^0\rightarrow K^-\pi^+\pi^0\pi^0$ decay. We investigate the sub-structures and determine the relative fractions and the phases among the different intermediate processes. Our results are used to provide an accurate detection efficiency and allow measurement of ${\cal B}(D^0\rightarrow K^-\pi^+\pi^0\pi^0) \,=\, (8.86 \pm 0.13(\text{stat}) \pm 0.19(\text{syst}))\%$.


I. INTRODUCTION
Many measurements of D meson decays have been performed since the D mesons were discovered in 1976 by Mark I [1,2].Today, most of the low-multiplicity D decay mode branching fractions (BFs) are well-measured.The largest decay modes are Cabibbo-favored (CF) hadronic (semileptonic) decay modes resulting from c → sW + , W + → u d (l + ν l ) transitions, but some of these decays are still unmeasured, in which the D 0 → K − π + π 0 π 0 decay should be the largest unmeasured mode.Chargeconjugate states are implied throughout this paper.
The D 0 /D + meson is the lightest meson containing a single charm quark.No strong decays are allowed, which makes the D 0 /D + meson a perfect place to study the weak decay of the charm quark.The CF Kπ, K2π, and K3π modes are the most common hadronic decay modes of D 0 /D + mesons.All Kπ and K2π branching fractions have been measured, but only four of the seven K3π [3] have been determined.The Mark III and The E691 collaboration performed amplitude analyses of all four D → Kπππ decay modes, K − π + π + π − , K 0 S π + π + π − , K − π + π + π 0 , and K 0 S π + π − π 0 [4,5].Recently, BESIII has remeasured the structure of the D 0 → K − π + π + π − decay with better precision [6].However, K3π modes with two or more π 0 's remain unmeasured.
Furthermore, the D 0 → K − π + π 0 π 0 decay has a large BF and is often used as a D 0 meson "tag mode" in experiment, such as in the CLEO and BESIII studies of D 0 semileptonic decays [7,8].This mode contributes up to 10% of the total reconstructed tags.Therefore, the accurate measurement of its sub-structures and branching fraction is essential to reduce systematic uncertainties of such analyses.While it is true that tag-mode BFs and sub-structure effects cancel to first order, higher-order systematic effects are increasingly important as statistics and precision increase.
The BESIII detector collected a 2.93 fb −1 dataset in 2010 and 2011 at √ s = 3.773 GeV [9,10], which corresponds to the mass of the ψ(3770) resonance.The ψ(3770) decays predominantly to D 0 D0 or D + D − without any additional hadrons.The excellent tracking, precision calorimetry, and a large D D threshold data sample at BESIII provide an excellent opportunity for study of the unmeasured D 0 → K − π + π 0 π 0 decay mode.The knowledge of intermediate structure will be crucial for determining the detection efficiency and useful for future usage as a tagging mode.We report here the first partial wave analysis (PWA) and BF measurement of the D 0 → K − π + π 0 π 0 decay.

II. DETECTION AND DATA SETS
The BESIII detector is described in detail in Ref. [11].The geometrical acceptance of the BESIII detector is 93% of the full solid angle.Starting from the interaction point (IP), it consists of a main drift chamber (MDC), a time-of-flight (TOF) system, a CsI(Tl) electromagnetic calorimeter (EMC) and a muon system (MUC) with layers of resistive plate chambers (RPC) in the iron return yoke of a 1.0 T superconducting solenoid.The momentum resolution for charged tracks in the MDC is 0.5% at a transverse momentum of 1 GeV/c.Monte Carlo (MC) simulations of the BESIII detector are based on geant4 [12].The production of ψ(3770) is simulated with the kkmc [13] package, taking into account the beam energy spread and the initial-state radiation (ISR).The photos [14] package is used to simulate the final-state radiation of charged particles.The evtgen [15] package is used to simulate the known decay modes with BFs taken from the Particle Data Group (PDG) [16], and the remaining unknown decays are generated with the LundCharm model [17].The MC sample referred to as generic MC, including the processes of ψ(3770) decays to D D, non-D D, ISR production of low mass charmonium states and continuum (e + e − → e + e − , µ + µ − , γγ and q q) processes, is used to study the background contribution.The effective luminosities of the generic MC sample correspond to at least 5 times the data luminosity.The signal MC sample includes D 0 → K − π + π 0 π 0 versus D0 → K + π − events generated according to the results of the fit to data.

III. EVENT SELECTION
Photons are reconstructed as energy clusters in the EMC.The shower time is required be less than 700 ns from the event start time in order to suppress fake photons due to electronic noise or e + e − beam background.Photon candidates within | cos θ| < 0.80 (barrel) are required to have larger than 25 MeV energy deposition and those with 0.86 < | cos θ| < 0.92 (endcap) must have larger than 50 MeV energy deposition.To suppress noise from hadronic shower splitoffs, the calorimeter positions of photon candidates must be at least 10°away from all charged tracks.Charged track candidates from the MDC must satisfy | cos θ| < 0.93, where θ is the polar angle with respect to the direction of the positron beam.The closest approach to the interaction point is required to be less than 10 cm in the beam direction and less than 1 cm in the plane perpendicular to the beam.
Charged tracks are identified as pions or kaons with particle identification (PID), which is implemented by combining the information of dE/dx in the MDC and the time-of-flight from the TOF system.For charged kaon candidates, the probability of the kaon hypothesis is required to be larger than that for a pion.For charged pion candidates, the probability for the pion hypothesis is required to be larger than that for a kaon.
The π 0 candidates are reconstructed through π 0 → γγ decays, with at least one barrel photon.The diphoton invariant mass is required to be in the range of 0.115 < M γγ < 0.150 GeV/c 2 .Two variables, beam constrained mass M BC and energy difference ∆E, are used to identify the D meson, where | p D | 2 and E D are the total reconstructed momentum and energy of the D candidate in the center-of-mass frame of the ψ(3770), respectively, and E beam is the calibrated beam energy.The D signals will be consistent with the nominal D mass in M BC and with zero in ∆E.
After charged kaons and charged pions are identified, and neutral pions are reconstructed, hadronic D decays can be reconstructed with a DTag technique.There are two types of samples used in the DTag technique: single tag (ST) and double tag (DT) samples.In the ST sample, only one D or D meson is reconstructed through a chosen hadronic decay without any requirement on the remaining measured tracks and showers.For multiple ST candidates, only the candidate with the smallest |∆E| is kept.In the DT sample, both D and D are reconstructed, where the meson reconstructed through the hadronic decay of interest is called the "signal side", and the other meson is called the "tag side".For multiple DT candidates, only the candidate with the smallest summation of |∆E|s in the signal side and the tag side is kept.
In this amplitude analysis, the DT candidates used are required to have the D 0 meson decaying to K − π + π 0 π 0 as the signal, and the D0 meson decaying to K + π − as the tag.For charged tracks of the signal side, a vertex fit is performed and the χ 2 must be less than 100.To improve the resolution and ensure that all events fall within the phase-space boundary, we perform a three-constraint kinematic fit in which the invariant masses of the signal D candidate and the two π 0 's are constrained to their PDG values [16].The events with kinematic fit χ 2 > 80 are discarded.
The tag side is required to satisfy 1.8575 < M BC < 1.8775 GeV/c 2 and −0.03 < ∆E < 0.02 GeV.The signal side is required to satisfy 1.8600 < M BC < 1.8730 GeV/c 2 and −0.04 < ∆E < 0.02 GeV.A K 0 S → π 0 π 0 mass veto, M π 0 π 0 / ∈ (0.458, 0.520) GeV/c 2 , is also applied on the signal side to remove the dominant peaking background, D 0 → K − K 0 S π + .The M BC and ∆E distributions of the data and generic MC samples are given in Fig. 1, where the generic MC sample is normalized to the size of data.Note that we always apply the ∆E requirements before plotting M BC , and vice-versa.
The generic MC sample is used to estimate the background of the DT candidates in the amplitude analysis.The dominant peaking background arises from D 0 → K − K 0 S π 0 , which is suppressed by the K 0 S mass veto from 2.2% to 0.07%.The remaining non-peaking background is about 1.0%.With all selection criteria applied, 5,950 candidate events are obtained with a purity of 98.9%.

IV. AMPLITUDE ANALYSIS
This analysis aims to determine the intermediate-state composition of the D 0 → K − π + π 0 π 0 decay.This fourbody decay spans a five-dimensional space.The daughter particle momenta are used as inputs to the probability density function (PDF) which describes the distribution of signal events.This is then used in an unbinned maximum likelihood fit to determinate the intermediate-state

A. Likelihood Function Construction
The PDF is used to construct the likelihood of the amplitude mode and it is given by where p is the set of the four daughter particles' four momenta and a is the set of the complex coefficients for amplitude modes.The (p) is the efficiency parameterized in terms of the daughter particles' four momenta.The four-body phase-space, R 4 , is defined as where α indicates the four daughter particles.This analysis uses an isobar model formulation, where the signal decay amplitude, A(a, p), is represented as a coherent sum of a number of two-body amplitude modes: where a i is written in the polar form as ρ i e iφi (ρ i is the magnitude and φ i is the phase), and A i (p) is the amplitude for the i th amplitude mode modeled as where the indexes 1 and 2 correspond to the two intermediate resonances.Here, F D i (p) is the Blatt-Weisskopf barrier factor for the D meson, while P 1,2 i (p) and F 1,2 i (p) are propagators and Blatt-Weisskopf barrier factors, respectively.The spin factor S i (p) is constructed with the covariant tensor formalism [18].Finally, the likelihood is defined as where k sums over the selected events and N s is the number of candidate events.Consequently, the log likelihood is given by ln L =  Since the second term of Eq. ( 7) is independent of a and the normalization integration in the denominator of the first term can be approximated by a phase-space MC integration, one can execute an amplitude analysis without knowing efficiency in advance.The phase-space MC integration is obtained by summing over a phase-space MC sample, where N g,ph is the number of generated phase-space events and N s,ph is the number of selected phase-space events.This holds since the generated sample is uniform in phase space, while the nonuniform distribution after selection reflects the efficiency.
For signal MC samples, the amplitude squared for each event should be normalized by the PDF which generates the sample.The normalization integration using signal MC samples is given by where N MC is the number of the signal MC sample and a gen is the set of the parameters used to generate the signal MC sample, which is obtained from the preliminary results using the phase-space MC integration.We allow for possible biases caused by tracking, PID, and π 0 data versus MC sample efficiency differences by introducing the correction factors γ , where j,data and j,MC are the π 0 reconstruction, the PID, or the tracking efficiencies as a function of p for the data and the MC sample, respectively.By weighting each signal MC event with γ , the MC integration is given by (11)

Spin Factor
For a decay process of the form a → bc, we use p a , p b , p c to denote the momenta of the particles a, b, c, respectively, and r a = p b − p c .The spin projection operators [18] are defined as µµ (a)P νν (a) + P µν (a)P νµ (a)) µ ν (a) . ( The covariant tensors are given by t( 1) We list the 10 kinds of spin factors used in this analysis in Table I, where scalar, pseudo-scalar, vector, axial-vector, and tensor states are denoted by S, P , V , A, and T , respectively.

Blatt-Weisskopf Barrier Factors
The Blatt-Weisskopf barrier F i (p j ) is a barrier function for a two-body decay process, a → bc.The Blatt-Weisskopf barrier depends on angular momenta and the magnitudes of the momenta of daughter particles in the rest system of the mother particle.The definition is given by where L denotes the angular momenta, and z = qR with q the magnitudes of the momenta of daughter particles in the rest system of the mother particle and R the effective radius of the barrier.For a process a → bc, we define while the values of R used in this analysis, 3.0 GeV −1 and 5.0 GeV   [18].Scalar, pseudo-scalar, vector, axial-vector, and tensor states are denoted by S, P , V , A, and T , respectively.

Decay chain S(p) D[S] → V1V2
t(1)µ (V1) t( 1) on evtgen).However, these values will also be varied as a source of systematic uncertainties.The X L (q) are given by

Propagators
We use the relativistic Breit-Wigner function as the propagator for the resonances K * 0 , K * − , and a 1 (1260) + , and fix their widths and masses to their PDG values [16].The relativistic Breit-Wigner function is given by where m = E 2 − p 2 and m 0 is the rest mass of the resonance.Γ(m) is given by where q 0 indicates the value of q when s a = m 2 0 .Resonances K1 (1270) 0 and K 1 (1270) − are also parameterized by the relativistic Breit-Wigner function but with constant width Γ(m) = Γ 0 since these two resonances are very close to the threshold of ρK and Γ(m) vary very rapidly as m changes.We parameterize the ρ with the Gounaris-Sakurai lineshape [19], which is given by The function f (m) is given by where and The normalization condition at P GS (0) fixes the param-

Kπ S-Wave
The kinematic modifications associated with the Kπ S-wave are modeled by a parameterization from scattering data [20,21], which are described by a K * 0 Breit-Wigner along with an effective range non-resonant component with a phase shift, with where a and r are the scattering length and effective interaction length, respectively.The parameters F (φ F ) and R(φ R ) are the magnitude (phase) for non-resonant state and resonance terms, respectively.The parameters M , F , φ F , R, φ R , a, r are fixed to the results of the D 0 → K 0 S π + π − analysis by BABAR [21], given in Table II.Note that we have also tested different parametrizations of the ππ S-wave, but no significant improvement is observed.We decide to use phase-space for the ππ S-wave.

B. Fit Fraction
The fit fraction (FF) is independent of the normalization and phase conventions in the amplitude formalism, and hence provides a more meaningful summary of amplitude strengths than the raw amplitudes, ρ i in Eq. ( 4), alone.The definition of the FF for the i th amplitude is where the integration is approximated by a MC integration with a phase-space MC sample.Since the FF does not involve efficiency, the MC sample used here is at the generator level instead of at the reconstruction level, as shown previously in Eq. (8).
As for the statistical uncertainty of the FF, it is not practical to analytically propagate the uncertainties of the FFs from that of the amplitudes and phases.Instead, we randomly perturb the variables determined in our fit (by a Gaussian-distributed amount controlled by the fit uncertainty and the covariance matrix) and calculate the FFs to determine the statistical uncertainties.We fit the distribution of each FF with a Gaussian function and the width is reported as the uncertainty of the FF.

C. Results of Amplitude Analysis
We perform an unbinned likelihood fit using the likelihood described in Section IV A, where only the complex a i are floating.Starting with amplitude modes with significant contribution, we add (remove) amplitude modes into (from) the fit one by one based on their statistical significances, which are obtained by the change of the log-likelihood value ∆ ln L with or without the amplitude mode under study.There are 26 amplitudes each with a significance larger than 4σ chosen as the optimal set, listed in Table III and the uncertainties are discussed in Section VI A. There are more than 40 amplitudes tested but not used in the optimal set (< 4σ significance), listed in Appendix A.
The amplitude D → K − a 1 (1260) + , a 1 (1260) + → ρ + π 0 [S] is expected to have the largest FF.Thus, we choose this amplitude as the reference (phase is fixed to 0) in the PWA.Other important amplitudes are D → (K − π 0 ) S ρ + , D → K − a 1 (1260) + with a 1 (1260) + [S] → ρ + π 0 , and D → K − a 1 (1260) + with a 1 (1260) + [S] → ρ + π 0 .The notation [S] denotes a relative S-wave between daughters in a decay, and similarly for [P ], [D].A MC sample is generated based on the PWA results, called the PWA signal MC sample.The projections of the data sample and the PWA signal MC sample on the invariant masses squared and the cosines of helicity angles for the K − π + , K − π 0 , π + π 0 and π 0 π 0 systems are shown in Fig. 2. The helicity angle θ ij (i or j is K − , π + and π 0 ) is defined as the angle between the momentum vector of the particle i in the ij rest frame and the direction of the ij system in the D rest frame.There are clear K * (892) 0 and K * (892) − resonances around 0.796 GeV 2 /c 4 in the M 2 K − π + and M 2 K − π 0 projections, respectively, and a ρ + (770) resonance around 0.593 GeV 2 /c 4 in the M 2 π + π 0 projection.The gap in the M 2 π 0 π 0 projection is due to the K 0 S mass veto.A more detailed goodness-of-fit study is presented in the next section.The PWA signal MC sample improves the accuracy of the DT efficiency (needed to determine the BF), which is discussed in more detail in Section V C.

D. Goodness-of-Fit
While the one-dimensional projections of the data sample and the PWA signal MC sample shown in Fig. 2 look quite good, much information is lost in projecting down from the full five-dimensional phase space.It is thus desirable to have a more rigorous test of the fit quality.We have programmed a "mixed-sample method" for determining the goodness of our unbinned likelihood fit [22].According to the method, we can calculate the "T" value of the mixing of two samples, the expectation mean, µ T , and the variance, σ 2 T .From these values, we can calculate a "pull", (T − µ T )/σ T , which should distribute as a normal Gaussian function due to statistical fluctuations.The pull is expected to center at zero if the two samples come from the same parent PDF, and be biased toward larger values otherwise.In the case of our PWA fit, the pull is expected to be a little larger than zero because some amplitudes with small significance are dropped.In other words, adding more amplitudes into the model is expected to decrease the pull.
To check the goodness-of-fit of our PWA results, we calculate the pull of the T value of the mixing of the data sample and the PWA signal MC sample, and it is determined to be 0.97, which indicates good fit quality.

V. BRANCHING FRACTION
We determinate the BF of D 0 → K − π + π 0 π 0 using the efficiency based on the results of our amplitude analysis.

A. Tagging Technique and Branching Fraction
To extract the absolute BF of the D 0 → K − π + π 0 π 0 decay, we obtain the ST sample by reconstructing the D0 meson through the D0 → K + π − decay, and the DT sample by fully reconstructing both D 0 and D0 through the D 0 → K − π + π 0 π 0 decay and the D0 → K + π − decay as the signal side and the tag side, respectively.The ST yield is given by and the DT yield is given by where N D 0 D0 is the total number of produced D 0 D0 pairs, B tag(sig) is the BF of the tag (signal) side, and ε are the corresponding efficiencies.
The BF of the signal side is determined by isolating B sig such that

B. Fitting Model
The ST yield, N ST tag , is obtained by a maximumlikelihood fit to the M BC (K + π − ) distribution.A Crystal Ball (CB) function [25], along with a Gaussian, is used to model the signal while an ARGUS function [26] is used to model the background.The signal shape is where f is a fraction ranging from 0 to 1, µ G and σ G are the mean and the width of the Gaussian function, respectively.The CB function has a Gaussian core transitioning to a power-law tail at a certain point, and is given by where N is the normalization and α controls the start of the tail.The beam energy (end point of the ARGUS function) is fixed to be 1.8865 GeV, while all other parameters are floating.
The DT yield, N DT tag,sig , is obtained by a maximumlikelihood fit to the two dimensional (2-D) M BC (K − π + π 0 π 0 ) versus M BC (K + π − ) distribution for the signal and the tag side with a 2-D fitting technique introduced by CLEO [23].This technique analytically models the signal peak, and considers ISR and mispartition (i.e., where one or more daughter particles are associated with the incorrect D 0 or D0 parent) effects, which are nonfactorizable in the 2-D plane.In this fitting, the mass of ψ( 3770) is fixed to be 3.773 GeV and the beam energy is fixed to be 1.8865 GeV.

C. Efficiency and Data Yields
An updated MC sample based on our PWA results, called the PWA MC sample, is used to determine the efficiency.The PWA MC sample is the generic MC sample with the K − π + π 0 π 0 versus K + π − events replaced by the PWA signal MC sample.All event selection criteria mentioned in Section III are applied except the M BC requirements.The projections to the signal and the tag side of the fit to the M BC distributions of the DT of data are shown in Figs.3(a) and (b), respectively.The background peak in the projection to the signal (tag) side axis is caused by events with a correct signal (tag) and a fake tag (signal).The fit to the M BC distribution of the ST of data is shown in Fig. 3(c), where both the mean values of the Gaussian function and the CB function agree well with our expectation for the D 0 mass.The ST and DT data yields are determined to be 534, 581 ± 769 and 6, 101 ± 83, respectively.The ST and DT efficiencies based on the PWA MC sample are (66.01 ± 0.03)% and (8.39 ± 0.04)%, respectively.

VI. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties of the PWA and BF measurement are discussed in Sections VI A and VI B, respectively.

A. Uncertainties for Amplitude Analysis
The systematic uncertainties for our amplitude analysis are studied in four categories: amplitude model, background, experimental effects, and fit bias.The contributions from the different categories to the systematic uncertainties for the FFs and phases are given in Tables IV and V, respectively.The uncertainties of these categories are added in quadrature to obtain the total systematic uncertainties.The effects of the amplitude model arise from three possible sources: the Kπ S-wave model, the effective barrier radii, and the masses and widths of intermediate particles.To determine the systematic uncertainties due to the Kπ S-wave model, the fixed parameters of the model are varied according to the BABAR measurement uncertainties [20,21], listed in Table II.The effective barrier radius R is varied from 1.5 to 4.5 GeV −1 for intermediate resonances, and from 3.0 to 7.0 GeV −1 for the D 0 .The masses and widths of intermediate particles are perturbed according to their published uncertainties in the PDG.The consequent changes of fitting results are considered as the systematic uncertainties inherent in the amplitude model.
The effects of background estimation are separated into non-peaking background, and peaking background.The uncertainties associated with non-peaking background are studied by widening the M BC and ∆E windows on the signal side to increase the non-peaking background.The peaking background can be mostly removed by the K 0 S mass veto.However, this veto is also a source of uncertainties.Its uncertainty is studied by widening this veto from the nominal M π 0 π 0 / ∈ (0.458, 0.520) GeV/c 2 to M π 0 π 0 / ∈ (0.418, 0.542) GeV/c 2 .
The experimental effects are related to the acceptance difference between data and MC sample caused by π 0 reconstruction, tracking, and PID efficiencies, which weight the normalization of the signal PDF, Eq. (10).To estimate the uncertainties associated with the experimental effects, the amplitude fit is performed varying π 0 reconstruction, tracking and PID efficiencies according to their uncertainties, and the changes of the nominal results are taken as the systematic uncertainties.
The fit bias is tested with 200 pseudoexperiment samples generated based on the PWA model.The distribution of each FF or phase is fitted with a Gaussian function.The difference of the mean and the nominal value is considered as the uncertainty associated with fit bias.

B. Uncertainties for Branching Fraction
We examine the systematic uncertainties for the BF from the following sources: tag side efficiency, tracking, PID, and π 0 efficiencies for signal, K − π + π 0 π 0 decay (PWA) model, yield fits, K 0 S peaking background and the K 0 S mass veto, and doubly Cabibbo-suppressed (DCS) decay.
The efficiency for reconstructing the tag side, D0 → K + π − , should almost cancel, and any residual effects caused by the tag side are expected to be negligible.Unlike the case of the tag side, the reconstruction efficiency of the signal side does not cancel in the DT to ST ratio.This efficiency of the signal side is determined with the PWA signal MC sample.The mismatches of tracking, PID, and π 0 reconstruction between the data and MC samples, therefore, bring in systematic uncertainties.
One possible source of those uncertainties is that the momentum spectra simulated in the MC sample do not match those in data, if there are any variations in efficiency versus momentum.This effect, however, is expected to be small due to the PWA MC sample's successful modeling of the momentum spectra in data, as shown in Fig. 2. The major possible source of the π 0 reconstruction, tracking, and PID systematic uncertainties is that, although the momentum spectra in the MC sample and data follow each other well, the efficiency of the MC sample disagrees with that of data as a function of momentum.This disagreement results in taking a correctly weighted average of incorrect efficiencies.We have performed an efficiency correction and choose 0.6%, 0.5%, 0.3%, and 0.2% as the systematic uncertainties for the π 0 reconstruction, kaon tracking, pion tracking, and kaon/pion PID, respectively.The uncertainty of the π 0 reconstruction efficiency is investigated with the control sample of D 0 → K − π + π 0 decays and the uncertainties for charged tracks and PID are determined using the control sample of D + → K − π + π + decays, D 0 → K − π + TABLE IV.FF systematic uncertainties (in units of statistical standard deviations) for: (I) the amplitude model, (II) background, (III) experimental effects, and (IV) fit bias.The total uncertainty is obtained by adding all contributions in quadrature.To estimate the systematic uncertainty caused by the imperfections of the decay model, we compare our PWA model to another PWA model which only includes amplitudes with significance larger than 5σ.The relative shift on efficiency is less than 0.5%.We therefore assign 0.5% as the systematic for the effect of any remaining decay modeling imperfections on efficiency.
To get the uncertainty of the yield fit, we change the nominal ∆E window to a wider one, −0.05 < ∆E < 0.03 GeV, and the change of the BF is considered as the associated uncertainty.
The K 0 S mass veto can veto most K 0 S peaking background and reduce it to be only 0.07% of the total events.However, the peaking background simulation is not perfect and the K 0 S mass veto also removes about 13% of the signal events.Thus, we estimate the uncertainty by narrowing the veto from M π 0 π 0 / ∈ (0.458, 0.520) GeV/c 2 to M π 0 π 0 / ∈ (0.470, 0.510) GeV/c 2 , while the K 0 S peaking background increases from 0.07% to 0.15% and the BF change is 0.18% of itself.We take this full shift as the corresponding uncertainty.
The smooth ARGUS background level is about 1.0% in the signal region.In addition, the 2-D M BC (K − π + π 0 π 0 ) versus M BC (K + π − ) fit works well for the background determination.Thus, we believe the uncertainty of the background with such small size will be very small and neglect it.
Our tag and signal sides are required to have oppositesign kaons.This means that our DT decays are either both CF or both DCS.These contributions can interfere with each other, with amplitude ratios that are approximately known, but with a priori unknown phase.The fractional size of the interference term varies between approximately ±2|A DCS /A CF | 2 ±2 tan 4 θ C , where θ C is the Cabibbo angle (the square in the first term arises as one power from each decay mode in the cross-term).The TABLE V. Phase, φ, systematic uncertainties (in units of statistical standard deviations) for: (I) the amplitude model, (II) background, (III) experimental effects, and (IV) fit bias.The total uncertainty is obtained by adding all contributions in quadrature.
Systematic uncertainties on the BF are summarized in Table VI, where the total uncertainty is calculated by a quadrature sum of individual contributions.

ACKNOWLEDGMENTS
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

VIII. APPENDIX A: AMPLITUDES TESTED
The following is a list of all amplitude modes tested and found to have a significance smaller than 4σ.These are not included in the final fit set.

FIG. 1 .
FIG. 1.The (a) MBC and (b) ∆E distributions on the tag side.The (c) MBC and (d) ∆E distributions on the signal side.The (red) arrows indicate the requirements applied in the amplitude analysis.The (blue) solid lines indicate the MC sample, while the (black) dots with error bars indicate data.

FIG. 3 .
FIG.3.Fits to the MBC distributions of the DT of the data sample projected to the (a) signal side (K − π + π 0 π 0 ) and the (b) tag side (K + π − ) and fit to (c) the MBC distributions of the ST of the data sample, where the (black) dots with error bars are data, the (red) solid lines are the total fit, the (green) dashed lines are the signal, and (blue) dotted lines are the background.

TABLE I .
Spin factor for each decay chain.All operators, i.e. t, have the same definitions as Ref.

TABLE II .
[21]meters of Kπ S-wave, by BABAR[21], where the uncertainties include the statistical and systematic uncertainties.

TABLE VI .
D 0 → K − π + π 0 π 0 BF systematic uncertainties.The total uncertainty is obtained by adding all contributions in quadrature.China (NSFC) underContract No. 11835012; National Natural Science Foundation of China (NSFC) under Contracts Nos.11475185, 11625523, 11635010, 11735014; 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.U1532257, U1532258, 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; German Research Foundation DFG under Contract No. Collaborative Research Center CRC 1044; Istituto Nazionale di Fisica Nucleare, Italy; Koninklijke Nederlandse Akademie van Wetenschappen (KNAW) under Contract No. 530-4CDP03; Ministry of Development of Turkey under Contract No. DPT2006K-120470; National Science and Technology fund; The Knut and Alice Wallenberg Foundation (Sweden) under Contract No. 2016.0157;The Swedish Research Council; U. S. Department of Energy under Contracts Nos.DE-FG02-05ER41374, DE-SC-0010118, DE-SC-0012069; University of Groningen (RuG) and the Helmholtzzentrum fuer Schwerionenforschung GmbH (GSI), Darmstadt.