Measurement of $\cos{2\beta}$ in $B^{0} \to D^{(*)} h^{0}$ with $D \to K_{S}^{0} \pi^{+} \pi^{-}$ decays by a combined time-dependent Dalitz plot analysis of BaBar and Belle data

We report measurements of $\sin{2\beta}$ and $\cos{2\beta}$ from a time-dependent Dalitz plot analysis of $B^{0} \to D^{(*)} h^{0}$ with $D \to K_{S}^{0} \pi^{+} \pi^{-}$ decays, where the light unflavored and neutral hadron $h^{0}$ is a $\pi^{0}$, $\eta$, or $\omega$ meson. The analysis is performed with a combination of the final data sets of the \babar\ and Belle experiments containing $471 \times 10^{6}$ and $772 \times 10^{6}$ $B\bar{B}$ pairs collected at the $\Upsilon\left(4S\right)$ resonance at the asymmetric-energy B factories PEP-II at SLAC and KEKB at KEK, respectively. We measure $\sin{2\beta} = 0.80 \pm 0.14 \,(\rm{stat.}) \pm 0.06 \,(\rm{syst.}) \pm 0.03 \,(\rm{model})$ and $\cos{2\beta} = 0.91 \pm 0.22 \,(\rm{stat.}) \pm 0.09 \,(\rm{syst.}) \pm 0.07 \,(\rm{model})$. The result for the direct measurement of the angle is $\beta = \left( 22.5 \pm 4.4 \,(\rm{stat.}) \pm 1.2 \,(\rm{syst.}) \pm 0.6 \,(\rm{model}) \right)^{\circ}$. The last quoted uncertainties are due to the composition of the $D^{0} \to K_{S}^{0} \pi^{+} \pi^{-}$ decay amplitude model, which is newly established by a Dalitz plot amplitude analysis of a high-statistics $e^{+}e^{-} \to c\bar{c}$ data sample as part of this analysis. We find the first evidence for $\cos2\beta>0$ at the level of $3.7$ standard deviations. The measurement excludes the trigonometric multifold solution $\pi/2 - \beta = (68.1 \pm 0.7)^{\circ}$ at the level of $7.3$ standard deviations and therefore resolves an ambiguity in the determination of the apex of the CKM Unitarity Triangle. The hypothesis of $\beta = 0^{\circ}$ is ruled out at the level of $5.1$ standard deviations, and thus CP violation is observed in $B^{0} \to D^{(*)} h^{0}$ decays.

We report measurements of sin 2β and cos 2β from a time-dependent Dalitz plot analysis of B 0 → D ( * ) h 0 with D → K 0 S π + π − decays, where the light unflavored and neutral hadron h 0 is a π 0 , η, or ω meson. The analysis is performed with a combination of the final data sets of the BABAR and Belle experiments containing 471 × 10 6 and 772 × 10 6 BB pairs collected at the Υ (4S) resonance at the asymmetric-energy B factories PEP-II at SLAC and KEKB at KEK, respectively. We measure sin 2β = 0.80 ± 0.14 (stat.) ± 0.06 (syst.) ± 0.03 (model) and cos 2β = 0.91 ± 0.22 (stat.) ± 0.09 (syst.) ± 0.07 (model). The result for the direct measurement of the angle is β = (22.5 ± 4.4 (stat.) ± 1.2 (syst.) ± 0.6 (model)) • . The last quoted uncertainties are due to the composition of the D 0 → K 0 S π + π − decay amplitude model, which is newly established by a Dalitz plot amplitude analysis of a high-statistics e + e − → cc data sample as part of this analysis. We find the first evidence for cos 2β > 0 at the level of 3.7 standard deviations. The measurement excludes the trigonometric multifold solution π/2 − β = (68.1 ± 0.7) • at the level of 7.3 standard deviations and therefore resolves an ambiguity in the determination of the apex of the CKM Unitarity Triangle. The hypothesis of β = 0 • is ruled out at the level of 5.1 standard deviations, and thus CP violation is observed in B 0 → D ( * ) h 0 decays.

I. INTRODUCTION
Breaking of CP symmetry is a small physical effect with profound consequences. CP violation causes particles and antiparticles to behave differently [1][2][3]. Even if the effects are tiny, CP violation provides the only possibility to assign matter and antimatter in an absolute and convention-independent way [4]. As one of the Sakharov requirements [5] for baryogenesis, CP violation is a key ingredient to generate the asymmetry between matter and antimatter shortly after the big bang that governs our present matter-dominated universe. However, CP violation in the standard model (SM) of electroweak interactions is several orders of magnitudes too small to account for the observed baryon asymmetry of the universe [6,7]. This is a strong motivation to search for additional sources of CP violation in nature. In the SM, the origin of CP violation is the single irreducible complex phase in the three-family Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix [8,9]. Testing this prediction of the Kobayashi-Maskawa theory [9] was the main objective for the construction and operation of the first-generation asymmetric-energy B factory experiments BABAR at SLAC (USA) and Belle at KEK (Japan). BABAR and Belle discovered CP violation in the decays of neutral and charged B mesons [10][11][12][13] and experimentally confirmed the theory predictions in numerous independent measurements [14].
In particular, BABAR and Belle observed CP violation in the interference between the direct decays of neutral B mesons into CP eigenstates and the decays after B 0 -B 0 oscillations (referred to as "mixing-induced CP violation") for the "gold plated" decay mode 1 B 0 → J/ψK 0 S and other decays mediated byb →ccs transitions [15,16]. By performing time-dependent CP violation measurements ofb →ccs transitions, BABAR and Belle precisely determined the parameter sin 2β ≡ sin 2φ 1 . 2 The angle β of the CKM Unitarity Triangle is defined as arg [−V cd V * cb /V td V * tb ], where V ij denotes a CKM matrix element. The current world average measured from b →ccs transitions is sin 2β = 0.691 ± 0.017 [17], which corresponds to an uncertainty on the angle β of 0.7 • . However, inferring the CP -violating weak phase 2β from * Now at: University of Huddersfield, Huddersfield HD1 3DH, UK † Also at: Università di Sassari, I-07100 Sassari, Italy the measurements of sin 2β is associated with the trigonometric two-fold ambiguity, 2β and π − 2β (a four-fold ambiguity in β), and therefore to an ambiguity in the determination of the apex of the CKM Unitarity Triangle.
The trigonometric ambiguity can be resolved experimentally by the measurements of B meson decays that involve multibody final states. Decay modes such as B 0 → J/ψK 0 S π 0 [18,19], B 0 → D * + D * − K 0 S [20,21], B 0 → K 0 S K + K − [22,23], B 0 → K 0 S π + π − [24,25], and B 0 → D ( * ) h 0 with D → K 0 S π + π − decays (abbreviated as B 0 → K 0 S π + π − ( * ) D h 0 ) [26][27][28][29] enable measurements of cos 2β in addition to sin 2β. Although sin 2β is precisely measured, the experimental uncertainties on cos 2β are sizable. Currently, the most precise single measurement has an uncertainty of approximately ±0.36 on the value of cos 2β [29]. However, no previous single measurement has been sufficiently sensitive to establish the sign of cos 2β that would resolve the trigonometric ambiguity without any assumptions. The strongest constraint in the direct estimation of the angle β was obtained by a measurement of B 0 → K 0 S K + K − decays by BABAR [22], which could resolve the ambiguity at the level of 4.8 standard deviations. However, B 0 → K 0 S K + K − decays do not provide a theoretically clean probe for the CP -violating weak phase 2β and only provide access to an effective weak phase β eff , because at leading order B 0 → K 0 S K + K − decays are not mediated by tree-level amplitudes but by quantum-loop ("penguin") transitions.
An experimentally elegant and powerful approach to access cos 2β and to resolve the trigonometric ambiguity is provided by B 0 → D ( * ) h 0 with D → K 0 S π + π − decays [26][27][28][29], where h 0 ∈ {π 0 , η, ω} denotes a light unflavored and neutral hadron. The decay B 0 → D * ω is not considered in this analysis. As shown in Fig. 1, the B 0 → D ( * ) h 0 decay is mediated only by tree-level amplitudes, and to a good approximation only by colorsuppressed, CKM-favoredb →cud tree amplitudes. Additional contributions from color-suppressed and doubly Cabibbo-suppressedb →ūcd amplitudes carry different weak phases, but are suppressed by a factor of |V ub V * cd /V cb V * ud | ≈ 0.02 relative to the leading amplitudes, and can be neglected at the experimental sensitivity of the presented measurement. The D 0 → K 0 S π + π − decay exhibits complex interference structures that receive resonant and nonresonant contributions from a rich variety of intermediate CP eigenstates and quasi-flavorspecific decays to the three-body final state. If the variations of the relative strong phase as a function of the D 0 meson three-body Dalitz plot phase space are known for D 0 → K 0 S π + π − decays, then both sin 2β and cos 2β can be measured from the time evolution of the B 0 → K 0 S π + π − ( * ) D h 0 multibody final state [26]. In an e + e − → Υ (4S) → B 0 B 0 event, the timedependent decay rate of the B 0 → K 0 S π + π − ( * ) D h 0 signal decays depends on the D 0 and D 0 decay amplitudes as a function of the three-body Dalitz plot phase space and on the CP -violating weak phase 2β, and is proportional to: The symbol ∆t denotes the proper-time interval between the decays of the two B mesons produced in the Υ (4S) event. The factor q = +1 (−1) represents the b-flavor content when the accompanying B meson is tagged as a B 0 (B 0 ). The parameters τ B 0 and ∆m d are the neutral B meson lifetime and the mass difference between the physical eigenstates of neutral B mesons ("B 0 -B 0 oscillation frequency"), respectively. The quantity η h 0 = (−1, −1, +1) is the CP eigenvalue of h 0 = (π 0 , η, ω), and the variable L is the orbital angular momentum of the Dh 0 and D * h 0 system. The relation η h 0 (−1) L equals −1 for Dh 0 , and +1 for D * h 0 (h 0 = ω). In this analysis, we consider only D * → Dπ 0 decays, so an additional factor of −1 that should be included for D * → Dγ decays need not be considered [30]. The D 0 and D 0 decay amplitudes where the symbol p i represents the four-momentum of a final state particle i.
Eq. (1) assumes no CP violation in B 0 -B 0 mixing and no direct CP violation in B 0 → D ( * ) h 0 decays. In our previous time-dependent CP violation analysis combining BABAR and Belle data [31], we determined the parameter C that measures direct CP violation in two independent samples of B 0 → D ( * ) h 0 decays. Using D meson decays both to CP eigenstates D CP → K + K − , K 0 S π 0 , and K 0 S ω, and using the high-statistics control sample provided by the CKM-favored D 0 → K + π − decay mode, no evidence for direct CP violation was found in either case [31]. This justifies the assumption of no direct CP violation in B 0 → D ( * ) h 0 decays for the present measurement.
The last term in Eq. (1) can be rewritten as: Eq. (2) allows the measurement of sin 2β and cos 2β as independent parameters by a time-dependent Dalitz plot analysis of B 0 → D ( * ) h 0 with D → K 0 S π + π − decays. Although elegant and appealing, the measurements of sin 2β and cos 2β in B 0 → D ( * ) h 0 with D → K 0 S π + π − decays are experimentally challenging and technically demanding. The branching fractions of these B and D meson decays are low, at the O(10 −4 ) and O(10 −2 ) level, respectively. These decay modes have neutral particles in the final states that lead to large backgrounds and low reconstruction efficiencies. In addition, either a detailed D 0 → K 0 S π + π − decay amplitude model or other experimental knowledge of the relative strong phase as a function of the D 0 meson three-body Dalitz plot phase space is required as input to perform the time-dependent Dalitz plot analysis of Time-dependent Dalitz plot analyses of B 0 → D ( * ) h 0 with D → K 0 S π + π − decays have been previously performed separately by BABAR and Belle. However, neither experiment was sensitive enough to establish CP violation [27][28][29]. Some of the measurements obtained results far outside of the physical region of the parameter space [27], and used different D 0 → K 0 S π + π − decay amplitude models [27,28], which complicates the comparison or the combination of the individual results.
In this article, we present measurements of sin 2β and cos 2β by a time-dependent Dalitz plot analysis of B 0 → D ( * ) h 0 with D → K 0 S π + π − decays that combines the BABAR and Belle data samples, totaling 1.1 ab −1 collected at the Υ (4S) resonance. In a recent combined analysis of the related decay, B 0 → D ( * ) CP h 0 with D CP denoting neutral D mesons reconstructed as two-body CP eigenstates, we demonstrated the technical feasibility and the physical advantage of the simultaneous analysis of the data collected by the BABAR and Belle experiments [31]. In the present measurement, the benefit is two-fold: first, the combination of the BABAR and Belle data samples improves the achievable experimental precision by effectively doubling the statistics available for the measurement; second, the combined approach enables common assumptions and the same D 0 → K 0 S π + π − decay ampli-tude model to be applied simultaneously in the analysis of the data collected by both experiments. The approach of combining BABAR and Belle data enables unique experimental sensitivity beyond what would be possible by combining two independent measurements, in particular for cos 2β. We derive the D 0 → K 0 S π + π − decay amplitude model from the data by a Dalitz plot amplitude analysis of a high-statistics e + e − → cc data sample. This approach ensures full control over the construction and the propagation of uncertainties of the D 0 → K 0 S π + π − decay amplitude model, and thus enables further improvement of the experimental sensitivity and robustness of the measurement.
The approach of combining the existing data of the B factory experiments BABAR and Belle results in measurements from a data sample with an integrated luminosity of more than 1 ab −1 . Data samples of comparable size are otherwise only achievable by future heavy flavor experiments: for example, the next-generation, highluminosity B factory experiment Belle II [32], which is expected to collect a data sample of 1 ab −1 by the year 2020 and is designed to collect 50 ab −1 by 2025. As such, the approach of combining the data from the firstgeneration asymmetric-energy B factory experiments enables not only unique experimental precision, but also demonstrates the discovery potential of Belle II at an early phase of the experiment.
The paper is structured as follows: Sect. II introduces the BABAR and Belle detectors and discusses the data sets used in the present analysis. In Sect. III, the Dalitz plot amplitude analysis to determine the D 0 → K 0 S π + π − decay model from a high-statistics e + e − → cc data sample collected by Belle is described. Sect. IV presents the measurements of sin 2β and cos 2β by a time-dependent Dalitz plot analysis of B 0 → D ( * ) h 0 with D → K 0 S π + π − decays combining the BABAR and Belle data sets. In Sect. V, the significance of the obtained results is studied. Finally, Sect. VI concludes the paper. The paper is accompanied by a Letter in Physical Review Letters [33].

II. THE BABAR AND BELLE DETECTORS AND DATA SETS
The results presented in this paper are based on data collected with the BABAR detector at the PEP-II e + e − storage rings [34] operated at the SLAC National Accelerator Laboratory (Menlo Park, USA) and with the Belle detector at the KEKB e + e − storage rings [35] operated at the KEK High Energy Accelerator Research Organization (Tsukuba, Japan). At PEP-II, 3.1 GeV positrons collide on 9 GeV electrons, and at KEKB, 3.5 GeV positrons collide on 8 GeV electrons. The centerof-mass (c.m.) energy of both PEP-II and KEKB is 10.58 GeV, which corresponds to the mass of the Υ (4S) resonance. Due to the asymmetry of the beam energies, the Υ (4S) is produced with a Lorentz boost of βγ = 0.560 at BABAR and 0.425 at Belle, allowing measurement of the proper-time interval between the decays of the two B mesons produced in Υ (4S) decays from the displacement of their decay vertices. The design of BABAR and Belle as asymmetric-energy B factory experiments is crucial to enable time-dependent CP violation measurements of neutral B mesons, as in the analysis presented in this paper.
The BABAR and Belle detectors are large-solid-angle multipurpose magnetic spectrometers, and are described in detail elsewhere [36][37][38]. The BABAR detector consists of a five-layer, double-sided silicon vertex tracker (SVT), a 40-layer drift chamber (DCH), an internally reflecting ring-imaging Cherenkov detector (DIRC), and a CsI(Tl) crystal electromagnetic calorimeter (EMC) located within a super-conducting solenoid magnet that provides a 1.5 T magnetic field. The instrumented flux return (IFR) of the solenoid magnet consists of iron plates interleaved with resistive plate chambers and, in the later runs, limited streamer tubes to detect K 0 L mesons and to identify muons.
The Belle detector consists of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), a barrellike arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter comprised of CsI(Tl) crystals (ECL) located inside a super-conducting solenoid coil that provides a 1.5 T magnetic field. An iron flux return located outside of the coil is instrumented to detect K 0 L mesons and to identify muons (KLM). Two inner detector configurations were used. A 2.0 cm radius beampipe and a 3-layer silicon vertex detector were used for the first sample of 152 × 10 6 BB pairs, while a 1.5 cm radius beampipe, a 4-layer silicon detector, and a smallcell inner drift chamber were used to record the remaining 620 × 10 6 BB pairs [39].
The Monte Carlo event generators used at BABAR and Belle are based on EvtGen [40], JETSET [41], and Photos [42]. The BABAR detector Monte Carlo simulation is based on Geant4 [43], and the Belle detector Monte Carlo simulation is based on Geant3 [44].
The first part of the analysis, described in Sect. III, is based on a data sample of 924 fb −1 recorded at or near the Υ (4S) and Υ (5S) resonances with the Belle detector [36]. This data set provides a high-statistics sample of e + e − → cc events that is used to determine the D 0 → K 0 S π + π − decay amplitudes. The data set provided by Belle enables a D 0 → K 0 S π + π − yield that is about three orders of magnitude larger than for the corresponding B meson decay to be studied by the combined BABAR+Belle approach. Therefore, the first part of the analysis does not require the combined use of the BABAR and Belle data sets.
The second part of the analysis, described in Sect. IV, is based on data samples collected at the Υ (4S) resonance containing (471 ± 3) × 10 6 BB pairs recorded with the BABAR detector and (772 ± 11) × 10 6 BB pairs recorded with the Belle detector. The combined BABAR and Belle data set is used to perform the time-dependent Dalitz plot analysis of B 0 → D ( * ) h 0 with D → K 0 S π + π − decays.
III. DETERMINATION OF THE D 0 → K 0 S π + π − DECAY AMPLITUDES BY DALITZ PLOT AMPLITUDE ANALYSIS USING BELLE e + e − → cc DATA

A. Event reconstruction and selection
The D * + → D 0 π + s candidates are reconstructed from D 0 → K 0 S π + π − decays and a low momentum ("slow") charged pion π + s . The slow pion enables the identification of the production flavor of the neutral D meson, which cannot be inferred directly from the self-conjugate three-body final state. The positive (negative) charge of the π + s determines the flavor of the neutral D meson to be D 0 (D 0 ). Neutral kaons are reconstructed in the decay mode K 0 S → π + π − , with the invariant mass required to be within 15 MeV/c 2 of the nominal value [45]. Further standard requirements exploiting the displacement of the K 0 S decay vertex from the interaction point (IP) described in Ref. [46] are applied. For candidates reconstructed from Υ (4S) and Υ (5S) data, requirements of p * (D * + ) > 2.5 GeV/c and p * (D * + ) > 3.1 GeV/c are applied, respectively, to reject combinatorial background and contamination from B meson decays, where p * denotes the momentum in the e + e − c.m. frame. The decay vertex of D * + candidates is determined by estimating the D 0 meson production vertex from a kinematic fit. In the kinematic fit, the D 0 meson is constrained to originate from the e + e − interaction region. The momentum resolution of soft pions is improved by a kinematic fit in which the π + s is constrained to the determined D * + decay vertex.
The reconstructed charmed meson decays are characterized by two observables: the D 0 candidate mass, M D 0 , and the D * + − D 0 mass difference, ∆M . Events are selected by requiring 1.825 < M D 0 < 1.905 GeV/c 2 and 140 < ∆M < 150 MeV/c 2 . For the Dalitz plot fit, a narrower, signal-enhanced region is defined by requiring B. Estimation of the D 0 → K 0 S π + π − signal and background yields The signal and background yields are estimated by a two-dimensional unbinned maximum-likelihood (ML) fit to the ∆M and M D 0 distributions. In the fit, the shape of the D * + → D 0 π + s with D 0 → K 0 S π + π − signal decays is parameterized by the sum of four two-piece normal distributions for M D 0 and by the sum of a normal distribution, a Johnson's SU function [47], a two-piece normal distribution, and a threshold function of the form (∆M − M π + ) 1/2 + a(∆M − M π + ) 3/2 + b(∆M − M π + ) 5/2 for ∆M . The width of the reconstructed ∆M distribution depends on the D 0 candidate mass. The ∆M distribution tends to become broader if the reconstructed D 0 mass deviates from the M D 0 peak position. To account for this correlation, the ∆M distribution is constructed by a conditional probability density function (p.d.f.) that scales the ∆M width with a fourth-order polynomial function that has the deviation of the reconstructed M D 0 from the M D 0 peak position as argument. In the fit, the fractions and widths of the tail components relative to that of the core components are fixed to values estimated using MC simulations, and the fractions and widths of the core components are determined by the fit.
The following four separate categories are considered for the background: The first source of background arises from the combination of correctly reconstructed D 0 → K 0 S π + π − candidates with random tracks during reconstruction. This "random slow pion" background has the same M D 0 shape as the signal, but the ∆M shape follows a smooth phase space distribution that is parameterized by a threshold function.
The second background category is composed of real π + s from D * + → D 0 π + s decays that are combined with wrong D 0 candidates formed from random tracks or with misreconstructed real D 0 decays. The distribution of this "real slow pion" background is mainly flat in M D 0 and very broad in ∆M due to the reconstruction of wrong D 0 candidates, but receives a small contribution that peaks in ∆M but is broad in M D 0 due to misreconstructed real D 0 decays. The shape of the background for wrong D 0 candidates is parameterized by a first-order polynomial function and a threshold function in M D 0 and ∆M , respectively; that for misreconstructed real D 0 decays is parameterized by a Crystal Ball function [48] and a Johnson's SU function for M D 0 and ∆M , respectively.
The third background category contains background from D 0 decay modes that have the same final state as D 0 → K 0 S π + π − decays, for example, D 0 → π + π − π + π − and D 0 → K 0 S K 0 S decays. The D 0 → π + π − π + π − decays are effectively removed by the applied K 0 S selection, and D 0 → K 0 S K 0 S decays have a very small branching fraction of O(10 −4 ). This "D 0 → 4π" background is parameterized by two Gaussian functions for M D 0 and the sum of a Gaussian function and a Johnson's SU function for ∆M . The yield of this background relative to the signal is at the sub-percent level. The fraction of this background is fixed to the expectation value obtained from Monte Carlo (MC) simulations.
The fourth background category accounts for the remaining combinatorial background originating from random combinations of tracks. This "combinatorial background" is parameterized by a first-order polynomial function in M D 0 and a threshold function in ∆M .
In the two-dimensional fit of the ∆M and M D 0 distributions, a total yield of 1 217 300 ± 2 000 signal events is obtained. The signal purity is 94% in the signal-enhanced with D 0 → K 0 S π + π − decays reconstructed from Belle e + e − → cc data. Component Yield D * + → D 0 π + s with D 0 → K 0 S π + π − signal 1 217 300 ± 2 000 Background containing real D 0 and random slow pions 61 330 ± 1 280 Background containing real slow pions and wrong D 0 249 700 ± 10 000 Background from D 0 → 4π 3 400 (fixed) Combinatorial background 271 000 ± 9 000 region used to extract the D 0 → K 0 S π + π − decay amplitude parameters. The results of the fit are summarized in Table I. The ∆M and M D 0 data distributions and projections of the fit are shown in Fig. 3.

C. Dalitz plot amplitude analysis
The D 0 → K 0 S π + π − decay proceeds via a rich variety of intermediate resonant and nonresonant modes contributing to the three-body final state. The contributions exhibit complex interference phenomena that are observ-able as characteristic patterns in the three-body Dalitz plot phase space as shown in Fig. 4. A Dalitz plot amplitude analysis is performed to disentangle and quantify the individual contributions.

Dalitz plot amplitude model
To describe the resonant and nonresonant substructure and to parameterize the D 0 → K 0 S π + π − decay amplitudes, the isobar ansatz [49] is combined with the Kmatrix formalism [50] for the ππ S-wave and the LASS   Events / 50 keV/c 2 parametrization [51] for the Kπ S-wave. In this approach, the D 0 → K 0 S π + π − decay amplitudes can be written as:

FIG. 3. (color online). Data distributions of ∆M and
The ππ and Kπ contributions with non-zero angular momentum are parameterized in the isobar ansatz by a coherent sum of the contributing intermediate quasi-twobody amplitudes. In the coherent sum, the r th interme-diate quasi-two-body amplitude A r enters with magnitude a r and relative phase φ r . The symbol F 1 denotes the decay amplitude for the ππ S-wave contributions parameterized by the K-matrix approach, and the symbol  A Kπ L=0 denotes the amplitude for the Kπ S-wave contribution using the LASS parametrization.
a. Isobar ansatz. In the isobar ansatz, the quasitwo-body amplitude for a neutral D meson decaying via the r th intermediate resonance (h 1 h 2 ) r with spin L to the three-body final state h 1 h 2 h 3 can be written as where the terms are described below. The form factors F (L) D and F (L) r describe the production D → rh 3 and the decay r → h 1 h 2 of the resonance r and the daughters of the resonance, respectively. The form factors are parameterized by Blatt-Weisskopf barrier penetration factors [52] that account for spindependent effects and prevent the decay amplitudes from diverging for large momentum transfers. The factors depend on the momentum q (p) of the bachelor particle h 3 (one of the resonance's daughter particles h 1 or h 2 ) evaluated in the resonance rest frame, and q 0 (p 0 ) is the value of q (p) when the invariant mass equals the pole mass of the resonance. The Blatt-Weisskopf barrier penetration factors are defined as where z = (|q|d) 2 and z 0 = (|q 0 |d) 2 . The parameter d represents the meson radius or the impact parameter of the decay particles for the D meson d D and the resonances d r , respectively. In the present analysis, The Zemach formalism [53] allows to describe the angular components of the amplitudes in a spin-tensor approach. The Zemach tensor formalism is applied to express the angular correlations among the final state particles by the function Z L (Ω), where the symbol Ω represent the angular relations of the involved particles.
The propagator term T r describes the dynamics in the resonance decay. In the present analysis, the term is parameterized by a relativistic Breit-Wigner (BW) lineshape function defined as where m 0 denotes the pole mass of the resonance, and the mass-dependent width Γ is given by The isobar ansatz is applied to parameterize the P -and D-wave contributions to the D 0 → K 0 S π + π − decay.
In the nominal Dalitz plot amplitude model, the following intermediate quasi-two-body resonances are included: the Cabibbo- To reduce the complexity of the Dalitz plot amplitude analysis, the masses and widths are fixed to the world averages [45] for all resonances except for the K * (892) ± , whose values are measured in the fit.
b. K-matrix formalism. The isobar ansatz has limitations, for example, in the case of broad and overlapping resonances or for resonances located close to thresholds of additional decay channels [49]. An alternative approach is provided by the K-matrix formalism, which preserves unitarity by construction in the presence of overlapping resonances and coupled channels. The K-matrix formalism is particularly suitable to describe the J P C = 0 ++ scalar contributions to the complex S-wave dynamics occurring in the π + π − system of D 0 → K 0 S π + π − decays. The BABAR, Belle, and LHCb experiments previously employed the K-matrix approach in Dalitz plot amplitude analyses of D 0 → K 0 S π + π − decays to perform measurements of D 0 -D 0 oscillations [54,55] and measurements of the Unitarity Triangle angle γ [56] in B meson decays [57,58]. Following the previous measurements, the K-matrix formalism in the P -vector approximation [59] is applied to model the ππ S-wave contribution to the D 0 → K 0 S π + π − decay. In this parametrization, the decay amplitude F 1 entering in Eq. (3) as the contribution of the ππ S-wave is defined by the relation where the indices i and j denote the particular channels (1 = ππ, 2 = KK, 3 = ππππ, 4 = ηη, and 5 = ηη ) contributing to the scattering process. The production vector P parameterizes the initial production of states into the open channels, and the K-matrix describes the scattering process. In this analysis, only the π + π − final states are considered, and s is the square of the invariant mass of the π + π − system. The terms I and ρ are the identity matrix and the phase-space matrix, respectively. The K-matrix is defined as The parameters m α are the physical poles of the Kmatrix, while g α i are the coupling constants of the i-th channel to the pole α. The parameters f scatt ij and s scatt 0 describe the smooth part of the K-matrix that is slowly varying. The unit of the number 1 is in GeV/c 2 . The symbol f A0 is the so-called "Adler zero" factor, defined as: This factor suppresses the false kinematic singularity at s = 0 in the physical region close to the π + π − threshold [60].
The production vector P has the same pole structure as the K-matrix and is defined as: The β α are the complex production couplings, and the parameters f prod 1j and s prod 0 describe the production of the slowly-varying part of the K-matrix.
In the present analysis, the K-matrix parameters m α , , s A0 , and s A are fixed to the results of a global analysis of available ππ scattering data [57,61] as summarized in Table II. The complex production couplings β α and the production parameters f prod 1j are free parameters determined from the fit.
c. LASS parametrization. For the Kπ S-wave, an approach introduced by the LASS collaboration to describe K − π + scattering processes is applied [51]. The Cabibbo-favored K * 0 (1430) − and the doubly Cabibbosuppressed K * 0 (1430) + contributions are each described by the empirical LASS parametrization. The LASS parametrization is constructed from a BW term for the K * 0 (1430) and a nonresonant component that has an effective range and introduces a phase shift: where The parameters R (φ R ) and F (φ F ) are the amplitudes (phases) of the resonant and nonresonant components, respectively. The parameters a and r are the scattering length and the effective interaction length, and q represents the momentum of the spectator particle in the Kπ rest frame. The parameters M and Γ(M 2 Kπ ) are the mass and the mass-dependent width of the resonant term defined in Eq. (9), and the phases δ R and δ F depend on m 2 Kπ . According to Ref. [57], this parametrization is equivalent to a K-matrix approach that describes a rapid phase shift originating from the resonant term and a slowly rising phase shift originating from the nonresonant term. The mass and the width of the K * 0 (1430) ± and the LASS R, φ R , F , φ F , a, and r are free parameters measured in the fit. The LASS parameters are required to be the same for the Cabibbo-favored K * 0 (1430) − and the doubly Cabibbo-suppressed K * 0 (1430) + contributions.

Dalitz plot reconstruction efficiency correction
Experimental effects, for example from the detector acceptance, the reconstruction algorithms, or the event selection, can induce non-uniformities for the reconstruction efficiency as a function of the Dalitz plot phase space, To account for these effects in the Dalitz amplitude analysis, the efficiency variations are estimated using a high-statistics sample of MC events of inclusive e + e − → cc decays that contain the D * + → D 0 π + s , In the MC simulations, the D 0 → K 0 S π + π − decay is generated uniformly in the available D meson decay phase space to uniformly populate the Dalitz plot. The generated decays are passed to a GEANT3-based simulation with a specific Belle configuration to simulate the detector response.
The simulated detector response then undergoes the same reconstruction algorithms and event selection requirements as for the data. The generated MC sample contains 50 × 10 6 D * + → D 0 π + s , D 0 → K 0 S π + π − signal decays, approximately 50 times the signal size in data, which enables the construction of a detailed map of the reconstruction efficiency.
The efficiency map is constructed using an approach BABAR introduced in the search for the Z(4430) − state [62]. In this approach, the efficiency is expressed as a function of the square of the two-body invariant mass M 2 K 0 S π − and cos θ K 0 S . The variable cos θ K 0 S is computed by the normalized dot product between the K 0 S π − three-momentum vector measured in the D meson rest frame and the three-momentum vector of the K 0 S meson after a Lorentz transformation from the D meson rest frame to the K 0 S π − rest frame. This choice of variables naturally introduces a "rectangular Dalitz plot" that is insensitive to potential binning effects that may arise at the curved edges of the M 2 K 0 S π − and M 2 π + π − Dalitz phase space due to the finite MC sample statistics. In order to parameterize the reconstruction efficiency and to smooth statistical fluctuations, the efficiency map is constructed in the following way.
In the first step, the angular variations of the efficiency are estimated by expanding the cos θ K 0 S distributions by a linear combination of Legendre polynomials up to order L = 7: The mass-squared dependent coefficients c L are estimated by fitting the linear combination of Legendre polynomials to the cos θ K 0 S distributions in intervals of The chosen order for the polynomial functions has been found to be sufficient to describe the details of the efficiency variations and at the same time to be low enough to avoid overfitting any structures. The dependence on the chosen order of the expansion in linear combinations of Legendre polynomials is weak, and lower or higher choices than L = 7 yield consistent results.
The reconstruction efficiency is almost flat over large parts of the Dalitz plot phase space. The efficiency decreases slightly at larger values of M 2 K 0 S π − and drops close to the kinematic border. The two-dimensional binned distributions of the reconstruction efficiency and the obtained parameterized efficiency maps are shown as a function of M 2 K 0 S π − and M 2 π + π − , and of M 2  Reconstruction Efficiency

Likelihood function and procedure for the
The D 0 → K 0 S π + π − decay amplitude parameters are estimated by an unbinned ML fit to the Dalitz plot distributions of the flavor-tagged D 0 sample. The likelihood function accounting for the contributions of the signal and background is written as where the index i runs over the reconstructed D 0 → K 0 S π + π − candidates. The signal fraction f sig and the fraction of the random slow pion background f rnd are determined by the two-dimensional fit to the M D 0 and ∆M distributions. The functions p sig , p rnd , and p bkg are the p.d.f.s of the Dalitz plot distributions for the signal, the random slow pion background, and the remaining background, respectively. The signal p.d.f. is constructed from the efficiency-corrected Dalitz plot intensities, computed from the absolute square of the Eq. (3), and by normalizing to the available Dalitz plot phase space: The random slow pion background is composed of a mixture of real D 0 and D 0 mesons decaying to the K 0 S π + π − final state. During the reconstruction of D * + → D 0 π + s decays, these D mesons are combined with random slow pion candidates. If the slow pion has the incorrect charge, the c-flavor content of the neutral D meson will be misidentified and the wrong flavor will be assigned. Neglecting possible production or detection asymmetries, the naïve expectation of the probability to select a slow pion track with the wrong charge is p = 0.5. The decay amplitudes for D 0 and D 0 mesons are related by an exchange of the Dalitz plot variables, The f wtag quantifies the fraction of "wrong D meson flavor-tags" and is estimated directly from the data by a separate Dalitz plot fit to the 150 < ∆M < 155 MeV/c 2 sideband region that has an enhanced population from the random slow pion background and no signal. In this Dalitz plot fit to the data sideband, the fraction of wrong D meson flavor-tag is measured and the result is f wtag = 0.492±0.075, in agreement with the naïve expectation. In the subsequent Dalitz plot fit to the signal region, f wtag is fixed to the estimate obtained from the sideband.
The background p.d.f. p bkg is constructed from the parameterized background model described in Sect. III C 3. The background is composed of combinatorial background and additional contributions from processes containing real slow pions and wrong D 0 mesons.
Due to the high statistics of the Belle e + e − → cc data sample of more than 10 6 events, and the complexity of the D 0 → K 0 S π + π − decay amplitude model, maximizing the likelihood function and performing the Dalitz plot fit is computationally intensive, taking hours to days on a single CPU core of a recent Intel Xeon processor-based Linux workstation. A new software framework for Dalitz plot amplitude analyses has been developed to increase the performance of the fit and to realize the present analysis. Key features of the framework are the parallel computing algorithms for both the evaluation of the likelihood function defined in Eq. (18), and for the numeric integration of the p.d.f.s. The parallel computing algorithms are realized using OpenMP [63,64] and enable the Dalitz plot fits to make simultaneous use of multiple CPUs to significantly reduce the required run time. In the present analysis, a speed-up of approximately a factor of 40 has been achieved for the time needed to reach convergence of the fit by using 64 CPU cores.
The D 0 → K 0 S π + π − decay amplitude parameters are determined by maximizing Eq. (18) for the Dalitz plot distributions in the signal-enhanced region defined in Sect. III A. The amplitude magnitudes a r and phases φ r of the intermediate resonant states are free parameters in the fit, and measured relative to the K 0 S ρ(770) 0 amplitude. The K 0 S ρ(770) 0 amplitude is fixed to a K 0 S ρ(770) 0 = 1 and φ K 0 S ρ(770) 0 = 0 • and serves as a reference. 5. Results of the D 0 → K 0 S π + π − Dalitz plot amplitude analysis The results for the estimated D 0 → K 0 S π + π − decay amplitude model parameters are summarized in Table III. The data distributions are shown in Figs. 4 and 6, and projections of the fit are shown in Fig. 6. The fit reproduces the data distributions well over the full range of the Dalitz plot. The fit projections exhibit few deviations, for example, for the ρ(770) 0 -ω(782) interference region in the M 2 π + π − projection. These deviations are very small compared to the overall scale of agreement.
The quality of the fit is estimated by a two-dimensional χ 2 test. The Dalitz plot data distributions are binned into square intervals with an edge length of 0.01 GeV/c 2 and then compared to the fit function. A reduced χ 2 of 1.05 is obtained for 31 272 d.o.f. based on statistical uncertainties only, indicating a good quality of the fit compared to previous models of this decay [54,55,57,65,66]. The normalized residuals contributing to the χ 2 function vary approximately uniformly over the Dalitz plot phase space and do not exhibit any macroscopic deviations or structures.
To quantify the contributions of individual amplitudes, the fit fractions (F F s) are evaluated. The F F for the r th intermediate resonant or nonresonant contribution is defined as: The sum of the fit fractions does not necessarily equal unity due to possible constructive or destructive interference effects among the amplitudes. In the present Dalitz plot amplitude analysis, the total fit fraction is 101.6%. The D 0 → K 0 S π + π − decay is dominated by the D 0 → K * (892) − π + mode which has a fit fraction of 59.9%. The second largest contribution is D 0 → K 0 S ρ(770) 0 with a fit fraction of 20.4%, followed by the π + π − S-wave with 10.0%.
To test further the agreement of the Dalitz plot amplitude model with the data, we follow an approach employed by BABAR in Ref. [67]. The Dalitz plot data distributions along the mass-squared directions are weighted by Y 0 k (cos θ) = (2k + 1)/4π P k (cos θ), where P k is the Legendre polynomial function of k th -order, and com-pared to the expectation of the corresponding Legendre moment computed from the Dalitz plot amplitude model. For M 2 K 0 S π − and M 2 π + π − , the weighted data distributions and the Legendre moments up to the 3 rd -order are shown in Fig. 7. The chosen representation is sensitive to the local phase and interference structures of the contributing amplitudes, complementary to the mass-squared projections. Good agreement is observed between the data distributions and the Dalitz plot amplitude model.

Model variations and crosschecks
The Dalitz plot amplitude analysis of D 0 → K 0 S π + π − decays is validated by various crosschecks. Before choosing the nominal Dalitz plot amplitude model, various alternative parameterizations and model variations have been considered.
The addition of further resonances (for example, the K * (1680) + π − mode) does not improve the fit quality nor result in significant fit fractions for these resonances. When parameterizing the ρ(770) 0 resonance by the Gounaris-Sakurai lineshape function [68] instead of the BW lineshape, worse agreement with the data is observed for the ρ(770) 0 and the ρ(770) 0 -ω(782) interference region. The determination of more parameters in the Dalitz plot fit (for example, the mass and the width of the ρ(770) 0 , ω(782), or other resonances) does not significantly improve the fit quality. In the nominal model, these parameters are fixed to the world averages [45] in order to reduce the complexity of the Dalitz plot fit.
Instead of the K-matrix and the LASS parametrization to describe the π + π − and Kπ S-waves, a model based on a pure isobar approach has been considered. In the isobar model, the π + π − S-wave is modeled by the σ 1 , σ 2 , f 0 (980), and f 0 (1370) resonances, and the Kπ Swaves by the Cabibbo-favored K * 0 (1430) − and the doubly Cabibbo-suppressed K * (1410) + resonances parameterized by BW lineshapes. An additional term that is constant in phase space is added to account for nonresonant contributions. For the isobar model, a reduced χ 2 of 1.23 is obtained for 31287 d.o.f. A similar isobar model including the σ 2 resonance has been used before by Belle [65,69] and CDF [66] in Dalitz plot amplitude analyses of D 0 → K 0 S π + π − decays. However, since the physical nature is not firmly established for all these states, in particular for the σ 2 resonance, and less agreement with the data was observed for the isobar model, it is not chosen as the nominal model.
The CLEO experiment performed a modelindependent determination of the relative strong phase between D 0 and D 0 → K 0 S π + π − decays by exploiting the quantum correlation of D 0 D 0 pairs produced from ψ(3770) decays in e + e − annihilations [70].  III. Results for the amplitude magnitudes ar, phases φr, fit fractions, K-matrix parameters for the π + π − S-wave, LASS parameters for the Kπ S-wave, and K * (892) ± parameters determined by the D 0 → K 0 S π + π − Dalitz plot fit performed for D * + → D 0 π + s events reconstructed from Belle data.
observed, corresponding to a p-value of 0.46. The results also agree well with a previous BABAR model of the same decay [57] that has been applied by CLEO to optimize the binning for the model-independent measurement of the relative strong phase.
IV. TIME-DEPENDENT DALITZ PLOT ANALYSIS OF B 0 → D ( * ) h 0 WITH D → K 0 S π + π − DECAYS USING BABAR AND BELLE DATA

A. Event reconstruction and selection
The similar performance of the BABAR and Belle detectors allows the use of almost identical selection require-ments in the two data sets. The event reconstruction and applied selection requirements discussed below follow the strategy used for the previous combined BABAR+Belle analysis of B 0 → D ( * ) CP h 0 decays described in Ref. [31]. Charged pion candidates are formed from tracks that are reconstructed from detected hits inside the tracking detectors and meet criteria for charged particles [36,37]. Photons are reconstructed from energy deposits of electromagnetic showers detected in the electromagnetic calorimeters. The energy of a photon candidate is required to be at least 30 MeV. in the decay modes η → γγ and π + π − π 0 . The invariant mass is required to be within [−25, +20] MeV/c 2 and ±10 MeV/c 2 of the nominal η mass [45] for η → γγ and η → π + π − π 0 candidates, respectively. The ω mesons are reconstructed in the decay mode ω → π + π − π 0 . The invariant mass of an ω meson candidate is required to be within [−15, +10] MeV/c 2 of the nominal ω mass [45].
Neutral kaons are reconstructed in the decay mode K 0 S → π + π − . The invariant mass of a K 0 S meson candidate is required to be within ±15 MeV/c 2 of the nominal value [45]. Standard selection requirements exploiting the displacement of the K 0 S decay vertex from the e + e − interaction point (IP) described in Refs. [46,71] are applied.
Neutral D mesons are reconstructed in the decay mode D → K 0 S π + π − . The invariant mass of a D meson candi-date is required to be within ±15 MeV/c 2 of the nominal value [45]. Neutral D * mesons are reconstructed in the decay mode D * → Dπ 0 . To select D * mesons, the reconstructed mass difference of neutral D * and D meson candidates is required to be within ±2.5 MeV/c 2 of the nominal value [45]. Neutral B mesons are reconstructed by combining light unflavored and neutral hadron candidates, h 0 ∈ {π 0 , η, ω}, with D ( * ) candidates. The decay modes B 0 → Dπ 0 , Dη, Dω, D * π 0 , and D * η, where sufficient signal yields are reconstructed, are included in the analysis. Neutral B mesons are selected using three variables that are constructed from kinematic observables: the beam-energy-constrained mass M bc , the energy difference ∆E, and the neural network classifier C NNout .
The beam-energy-constrained mass is defined as: where E * beam is the energy of either beam provided by the e + e − collider, the variables p * D ( * ) and E * D ( * ) are the three-momentum and the energy of the D ( * ) meson candidates, and p * h 0 and M h 0 are the three-momentum and the invariant mass of the h 0 candidates. Observables marked with an asterisk are evaluated in the e + e − c.m. frame. Belle introduced the variable M bc in the measurements of B meson decays mediated by radiative penguin transitions [72] as an alternative to the more commonly used variable We note that M bc does not directly depend on the three-momentum magnitude nor the energy, but only on the direction of flight of the h 0 candidate. Therefore, M bc is insensitive to potential correlations with the energy difference, defined as In the present analysis, non-trivial correlations emerge between M bc and ∆E for final states containing photons from the reconstructed h 0 decay modes due to energy mismeasurements by the electromagnetic calorimeters, for example, caused by shower leakage effects. The use of M bc effectively eliminates these correlations and enables factorizing the p.d.f.s constructed from the M bc and ∆E observables in multi-dimensional fits. The neural network combines information characterizing the shape of the events and is based on 16 modified Fox-Wolfram moments [73,74]. Following an approach introduced by Belle in Ref. [75], the variable C NNout is constructed from the output of the neural network classifier, C NNout , by the following transformation: The variables C min NNout and C max NNout are adjustable parameters, and are related to the output domain of C NNout . In this analysis, C min NNout = 0.2 and C max NNout = 1 are chosen. After the transformation to C NNout , the output of the neural network classifier exhibits smooth distributions around a peak position that differs for e + e − → qq (q ∈ {u, d, s, c}) continuum events and BB events. Candidates from continuum events tend to be distributed around a peak position at negative values of C NNout , while BB events are distributed around a peak position at positive values. The C NNout distributions can be described by empirical parameterized models with few d.o.f., such as the the Novosibirsk function, an empirical p.d.f. inspired by the log-normal distribution and defined in Ref. [76]. The use of a parameterized model has technical advantages when including the neural network classifier in addition to M bc and ∆E in multi-dimensional fits to extract the B 0 → D ( * ) h 0 signal. Before applying the transformation described above, a loose requirement of C NNout > 0.2 is applied to remove regions that are almost exclusively populated by continuum background events.
The following requirements are applied on M bc , ∆E, and C NNout to select neutral B mesons: 5.24 < M bc < 5.29 GeV/c 2 , −150 < ∆E < 200 MeV, and −8 < C NNout < 10. For B 0 → D ( * ) h 0 signal decays, the M bc , ∆E, and C NNout distributions exhibit smooth peaking structures. The shapes of the signal component are parameterized by two Novosibirsk functions for M bc , one symmetric and two two-piece normal distributions for ∆E, and two Novosibirsk functions for C NNout . The signal shapes are calibrated using the high-statistics data control sample of B 0 → D ( * )0 h 0 decays with the CKM-favored For B 0 → Dh 0 decays, candidates can originate from the corresponding B 0 → D * h 0 decay modes, if the slow neutral pion from D * → Dπ 0 decays is missed during the reconstruction. This "crossfeed component" originates from true B 0 → D * h 0 signal decays and has therefore signal-like properties. The crossfeed has similar shapes as the signal but peaks at negative ∆E. The contribution of the crossfeed is small, at the level of 3−13% with respect to the signal. In the fits, the fractions of this component are fixed to the values estimated from high-statistics MC simulations of signal decays. The shapes of the crossfeed component are parameterized by two Novosibirsk functions for M bc , one kernel density estimator for ∆E, and two Novosibirsk functions for C NNout .
In addition to the contributions from the signal and the signal-like crossfeed, the fit model accounts for the following three separate sources of background. The first source originates from partially-reconstructed B + → D ( * )0 ρ + decays, which constitute a background for B 0 → D ( * ) π 0 decays when the charged pion from ρ + → π + π 0 decays is soft. This background arises only for B 0 → Dπ 0 and B 0 → D * π 0 decays, but is not present for the other B 0 → D ( * ) h 0 decay modes. Like the crossfeed component, the background from B + → D ( * )0 ρ + decays has a similar shape as the signal, but peaks at negative ∆E. The shapes are parameterized by two Novosibirsk functions for M bc , one kernel density estimator for ∆E, and two Novosibirsk functions for C NNout . The B + →D ( * )0 ρ + background component is determined by the fit.
The second source of background arises from B meson candidates formed from random combinations of final state particles originating from e + e − → BB events. This "combinatorial BB background" is low in the present analysis. The combinatorial BB background exhibits smooth phase space distributions in M bc and ∆E, and peaks at positive C NNout . The shapes are parameterized by an ARGUS function [77] for M bc , a second-order polynomial function for ∆E, and two Novosibirsk functions for C NNout .
The third source of background originates from e + e − → qq (q ∈ {u, d, s, c}) continuum events. This continuum background exhibits smooth phase space distributions in M bc and ∆E, and peaks at negative C NNout . The shapes are parameterized by an ARGUS function for M bc , a second-order polynomial function for ∆E, and two Novosibirsk functions for C NNout .
In total, B 0 → D ( * ) h 0 signal yields of 1 129 ± 48 events for BABAR and 1 567 ± 56 events for Belle are obtained. The signal yields separated by experiment and decay mode are summarized in Table IV. The experimental M bc , ∆E, and C NNout distributions and projections of the fits are shown in Fig. 8.

Decay mode BABAR
Belle

C. Time-dependent Dalitz plot analysis
At BABAR (Belle) the Υ (4S) is produced with a Lorentz boost of βγ = 0.560 (0.425), allowing the measurement of the proper-time interval between the decays of the two B mesons. The proper-time interval ∆t is given by ∆z/cβγ, where ∆z denotes the spatial distance between the decay vertices of the two B mesons in the laboratory frame. The BABAR and Belle techniques to measure the flavor-tagged proper-time intervals of the B mesons and to extract the CP violation parameters are described in detail in Refs. [14][15][16][78][79][80][81]. The B 0 → D ( * ) h 0 signal decay vertices are reconstructed by kinematic fits that include experimental knowledge of the IP position. For BABAR, the applied vertex reconstruction algorithm simultaneously includes the complete B meson decay tree, including all secondary decays, in the kinematic fit. For Belle, the vertex reconstruction is performed in an iterative bottom-up approach starting with the final state particles. The decay vertex and the b-flavor content of the accompanying B meson are estimated from the reconstructed decay products not assigned to the signal B meson. The b-flavor content is inferred by the flavortagging procedures described in Refs. [15,80]. The applied algorithms account for different signatures such as the presence and properties of prompt leptons, charged kaons, and pions originating from the decay of the accompanying B meson, and assign a flavor and an associated probability.
The experimental conditions and the instrumentation of the detectors are different for BABAR and Belle. The finite experimental resolution in the measurements of proper-time intervals are different for BABAR and Belle, and both experiments follow different approaches to describe the resolution effects. The two experiments employ different multivariate techniques for the flavor-tagging. BABAR uses a neural network-based approach and Belle uses a multi-dimensional likelihood approach.
The time-dependent Dalitz plot analysis to measure the CP violation parameters follows the technique established in the previous combined BABAR+Belle timedependent CP violation measurement of B 0 → D ( * ) CP h 0 decays [31]. The strategy of the combined approach is to apply established, experiment-specific techniques to describe proper-time resolution and flavor-tagging effects by BABAR and Belle to the data collected by the particular experiment. The combined measurement is then performed by maximizing the log-likelihood function constructed from the p.d.f.s and the data collected by both experiments: The indices i and j run over events reconstructed from BABAR and Belle data, respectively. All events used in the M bc , ∆E, and C NNout fits are included. The P are the p.d.f.s of the experimental flavor-tagged proper-time interval and Dalitz plot distributions of the B mesons measured in the events, and are defined as: The index k represents the signal and background components. The fractions of the components, f k , are evaluated on an event-by-event basis as a function of M bc , ∆E, and C NNout . The P k are the p.d.f.s that describe the particular underlying particle physics process and are the same for both experiments. The P k are convolved with the resolution functions R k that account for the finite proper-time resolution. For the signal, the p.d.f.s are constructed from Eqs. (1) and (2) convolved with the experiment-specific resolution functions to account for the finite proper-time resolution [15,79], and include the effect of incorrect flavor assignments by the applied flavor-tagging algorithms [15,80] and a correction to account for the variations of the reconstruction efficiency as a function of the position on the Dalitz plot.
Neutral D mesons produced in B 0 → D ( * ) h 0 decays have a different momentum spectrum than those produced in e + e − → cc events. In addition, the yield for the B 0 → D ( * ) h 0 decay modes studied by the combined BABAR+Belle approach is about three orders of magnitude lower than that for the D 0 → K 0 S π + π − decays reconstructed from e + e − → cc events. Therefore, the Dalitz plot reconstruction efficiency correction used for the analysis of B 0 → D ( * ) h 0 decays is different from that described in Sect. III C 2, and a parametrization with fewer d.o.f. is chosen. The reconstruction efficiency map is constructed separately for BABAR and Belle by the fit of a two-dimensional 3 rd -order polynomial function in the Dalitz plot variables M 2 K 0 S π − and M 2 K 0 S π + to the reconstruction efficiency distributions obtained from high-statistics samples of MC events of B 0 → D ( * ) h 0 with D → K 0 S π + π − signal decays. For the signal-like crossfeed from partiallyreconstructed B 0 → D * h 0 decays, the p.d.f.s are constructed as for the signal, but account for distinct properties such as the CP -eigenvalues of the particular final states of the crossfeed contribution. The charged B meson background from partially-reconstructed B + → D ( * )0 ρ + decays is parameterized by an exponential p.d.f. accounting for the B + lifetime convolved with the experiment-specific resolution functions. The combinatorial BB background and the continuum background share the same parametrization for BABAR and Belle. For each background component, the p.d.f.s are constructed from the sum of a Dirac delta function to model background from prompt particles and an exponential p.d.f. with effective lifetimes to model the non-prompt background. The background p.d.f.s are convolved with a resolution function modeled as the sum of two Gaussian functions whose widths depend linearly on the uncertainty of ∆t. The ∆t parameters for the combinatorial BB background and the continuum background are determined by fits to the M bc < 5.26 GeV/c 2 sidebands and are fixed in the measurement.
The linear correlation between sin 2β and cos 2β is 5.1%. The result deviates less than 1.0 standard deviations from the trigonometric constraint given by sin 2 2β + cos 2 2β = 1. An alternative fit is performed to measure directly the CP -violating phase β using the signal p.d.f. constructed from Eq. (1). The result of this fit is: The evaluation of the experimental systematic uncertainties and the uncertainties due to the applied D 0 → K 0 S π + π − decay amplitude model are described in Sects. IV D 1 and IV D 2.
The B 0 → K 0 S π + π − ( * ) D h 0 decays proceeds via a rich variety of intermediate CP eigenstates and quasiflavor-specific decays contributing to the multibody final state. These intermediate contributions involve different physics in the time evolution of the B meson decay, and hence exhibit different proper-time interval distributions. In Fig. 9, the proper-time interval distributions and projections of the fit for sin 2β and cos 2β are shown for two different regions of the D 0 → K 0 S π + π − phase space. In Figs. 9a and c, a region of phase space predominantly populated by CP eigenstates,  [82] prevents the two neutral B mesons from being produced with the same flavor at ∆t = 0, which in Figs. 9b and d is additionally smeared by experimental resolution effects. The time evolution follows a 1 ± cos(∆m∆t) distribution, and the corresponding time-dependent oscillation asymmetry exhibits a cosine oscillation.
Various cross-checks are performed to validate the procedure of the measurement. The B 0 → D ( * )0 h 0 decays with the CKM-favored D 0 → K + π − decay have very similar kinematics and background composition as B 0 → D ( * ) h 0 with D → K 0 S π + π − decays and provide a high-statistics control sample. In total, signal yields of 3 029±73 events for BABAR and 4 042±84 events for Belle are obtained for the control sample. Using the same analysis approach, the time-dependent CP violation measurement of the control sample yields both mixing-induced and direct CP violation consistent with zero, in agreement with the expectation of negligible CP violation for these flavor-specific decays. Measurements of the neutral B meson lifetime for B 0 → D ( * ) h 0 with D → K 0 S π + π − decays and for the control sample without flavortagging applied yield τ B 0 = (1.500 ± 0.052 (stat.)) ps and τ B 0 = (1.535 ± 0.028 (stat.)) ps, respectively, and are in agreement with the world average τ B 0 = (1.520 ± 0.004) ps [17]. In addition, all measurements have been performed for data separated by experiments and yield consistent results. The results for B 0 → D ( * ) h 0 with D → K 0 S π + π − decays separated by experiments are sin 2β = 0.91 ± 0.20 (stat.), cos 2β = 0.87 ± 0.31 (stat.), and β = (25.6 ± 6.4 (stat.)) • for BABAR, and sin 2β = 0.70 ± 0.20 (stat.), cos 2β = 0.96 ± 0.30 (stat.), and β = (19.6 ± 6.1 (stat.)) • for Belle, respectively.

D. Determination of the systematic uncertainties
The present analysis accounts for two classes of systematic uncertainties on the measured CP violation parameters: first, the experimental systematic uncertainty accounts for experimental effects that can affect the timedependent Dalitz plot analysis; second, the Dalitz plot model uncertainty accounts for assumptions made on the applied D 0 → K 0 S π + π − decay amplitude model.

Experimental systematic uncertainties
The estimation of the experimental systematic uncertainty on the CP violation parameters follows established methods, described in Refs. [15,16,31]. The evaluation of the individual contributions to the experimental systematic uncertainty are described below, and the results are summarized in Table V.
The systematic uncertainty due to vertex reconstruction accounts for the applied vertex reconstruction algorithms, the requirements applied to select B mesons, the uncertainty of the z scale, possible ∆t biases, and effects due to possible misalignment of the silicon vertex detectors. For the vertex reconstruction algorithms, the constraints in the kinematic fits and the applied selection requirements of the signal B meson and the accompanying B meson are varied. For BABAR, the uncertainty due to the z scale and the Lorentz boost is estimated by variations of the corresponding scale and uncertainties. For Belle, a possible ∆t bias is estimated using MC simulations. Misalignment effects of the silicon vertex detectors are studied by MC simulations, and corresponding systematic uncertainties are assigned.
Experiment-specific resolution models are applied to account for effects due to the finite experimental ∆t resolution. The ∆t resolution function parameters are fixed to values obtained from control samples using BABAR and Belle data. The systematic uncertainty due to the applied ∆t resolution functions is estimated by variation of the resolution model parameters within their uncertainties.
The parameters of the ∆t model for the combinatorial BB background and the continuum background are 0 10 20 Events / 1.5 ps (a) BABAR  Events / 1.5 ps (d) Belle  [83] with weights obtained from the fit presented in Fig. 8. Two different regions of the D → K 0 S π + π − phase space are shown. In the plots of the left column, a region predominantly populated by CP eigenstates, B 0 → K 0 S ρ(770) 0 ( * ) D h 0 , is selected by requiring |M ρ(770) − M π + π − | < 150 MeV/c 2 . In the plots of the right column, a region predominantly populated by quasi-flavor-specific determined by fits to the M bc < 5.26 GeV/c 2 data sidebands. The systematic uncertainty due to the background ∆t p.d.f.s is estimated by variation of the ∆t background model parameters within their uncertainties. The signal purity is estimated by the three-dimensional unbinned ML fit to the M bc , ∆E, and C NNout distributions. The systematic uncertainty due to the signal purity estimation is estimated by variation of the fit parameters within their uncertainties.
The b-flavor content of neutral B mesons is inferred by multivariate BABAR-and Belle-specific flavor-tagging algorithms. The flavor-tagging algorithms are calibrated using control samples reconstructed from BABAR and Belle data. The systematic uncertainty due to the flavortagging is estimated by variation of the wrong-tag fractions and the corresponding wrong-tag fraction differences for each tagging category within their uncertainties.
The neutral B lifetime τ B 0 , the charged B meson lifetime τ B + , and the B 0 -B 0 oscillation frequency ∆m d are fixed to the world averages. The systematic uncertainty due to these fixed physics parameters is estimated by variation of the lifetimes and oscillation frequency within their uncertainties. The systematic uncertainty due a possible small fit bias in ∆t measurements is estimated by MC simulations. Large MC samples are generated using a complex D 0 → K 0 S π + π − decay amplitude model and with CP violation, the same reconstruction algorithms and event selection requirements are applied to the MC samples as for the data, and the time-dependent Dalitz plot analysis is performed. The deviations of the central values of the CP violation parameters measured using the MC samples from the nominal result are assigned as systematic uncertainties.
The effect due to the applied Dalitz plot reconstruction efficiency correction for neutral D mesons produced in B 0 → D ( * ) h 0 decays is estimated by removing the efficiency correction. The time-dependent Dalitz plot analysis is performed without the efficiency correction, and assigning the deviations from the nominal result as systematic uncertainty due to the Dalitz plot reconstruction efficiency correction.
Most systematic uncertainties are independent for BABAR and Belle. Possible correlations such as for the physics parameters are considered. Additional contributions to the systematic uncertainty from possible sources of peaking background and the tag-side interference have been considered and can be neglected in the presented measurement.
The total experimental systematic uncertainty is the quadratic sum of all contributions.

Uncertainty due to the Dalitz plot amplitude model
The model uncertainty accounts for the dependence of the CP violation parameters on the D 0 → K 0 S π + π − de-cay amplitude model determined by the Dalitz plot amplitude analysis using the high-statistics Belle e + e − → cc data sample described in Sect. III C. The strategy to estimate the model uncertainty is to repeat the D 0 → K 0 S π + π − Dalitz plot amplitude analysis with alternative assumptions and variations of the D 0 → K 0 S π + π − decay amplitude model. The time-dependent Dalitz plot analysis of B 0 → D ( * ) h 0 decays is then performed using the alternative models as input, and the deviations from the result using the nominal D 0 → K 0 S π + π − decay amplitude model are assigned as model uncertainty on the CP violation parameters. The evaluation of the individual contributions to the uncertainty due to the Dalitz plot amplitude model are described below, and the results are summarized in Table VI.
For the masses and widths of resonances fixed to the world averages, each resonance parameter is varied within its uncertainty to estimate the associated model uncertainty.
The model uncertainty due to the chosen π + π − Swave parametrization using the K-matrix formalism is estimated by replacing the nominal K-matrix solution by alternative solutions from Ref. [61]. In addition, the parameter s prod 0 is varied within its uncertainty, which is taken from Ref. [57].
The LASS parametrization is used to model the Kπ S-waves. The model uncertainty is estimated by replacing the LASS parametrization for the K * 0 (1430) − and K * 0 (1430) + resonances by standard relativistic BW terms.
The model uncertainty due to the chosen Blatt-Weisskopf barrier factors for D mesons and intermediate resonances is estimated by varying the fixed parameters d D and d r each by ±0.5 c/GeV.
The fraction of wrong D meson flavor-tags of the flavor-tagged cc data sample is fixed to the value estimated from the fit to the ∆M sideband region on data. The D meson mistag fraction is varied within its uncertainty to evaluate the associated model uncertainty.
The model uncertainty due the applied Dalitz plot reconstruction efficiency correction is estimated by replacing the parameterized efficiency map by the corresponding two-dimensional binned distributions.
In the Dalitz plot amplitude analysis, the background is described by a parameterized model taken from the ∆M and M D 0 sideband regions on data. The model uncertainty due to the applied background description is estimated by replacing the parameterized background model by the two-dimensional binned distributions from the data sidebands.
Most intermediate two-body resonances contributing to D 0 → K 0 S π + π − decays have a natural width much larger than the finite experimental resolution of reconstructed invariant masses, and resolution effects can be neglected in the D → K 0 S π + π − Dalitz plot amplitude analysis. The ω(782) width, 8.5 MeV, is comparable to the mass resolution. To estimate the size of possible effects due to the mass resolution and to evaluate the asso- ciated model error, the width of the ω(782) is increased by 20%. The signal and background fractions used in the Dalitz plot amplitude analysis are determined by the fit of the two-dimensional ∆M and M D 0 distributions. The model uncertainty due to the signal purity estimation is determined by varying the the ∆M -M D 0 model parameters within their uncertainties.
The statistical uncertainties on the Dalitz plot amplitude model parameters that are summarized in Table III are caused by the finite size of the cc data sample. To propagate the statistical uncertainties to the CP violation parameters and assign the associated model error, each parameter is varied within its uncertainty. For individual resonances, the correlations between phases and amplitudes are accounted for. An explicit treatment of additional correlations between resonances important in the CP violation measurement were found to be negligible. The chosen approach has to be found sufficient given that this systematic uncertainty does not limit the precision of the measurement.
The dependence of the model on resonances with very small contributions is estimated by removing resonances with fit fractions of 0.1% or lower. The doubly Cabibbosuppressed K * (1410) + , K * 2 (1430) + , and K * 0 (1430) + , and the K * (1410) − are each removed from the model. For each model variation, the D 0 → K 0 S π + π − Dalitz plot amplitude analysis is repeated to estimate the associated model uncertainty.
As a further cross-check and estimate of the possible model-dependence, a pure isobar D 0 → K 0 S π + π − decay model without the K-matrix parametrization is constructed. As in the isobar model discussed in Sect. III C 6, the intermediate resonant contributions to the π + π − Swave are modeled by the σ 1 , σ 2 , f 0 (980), and f 0 (1370) resonances, and a term constant in phase space is included to account for nonresonant contributions. The D 0 → K 0 S π + π − Dalitz plot amplitude analysis and the time-dependent Dalitz plot analysis of B 0 → D ( * ) h 0 decays are repeated using the alternative model, and the deviations of the CP violation parameters from the baseline result are assigned as model uncertainty. The result with the isobar model agrees well with the baseline result, which indicates small overall model dependence and robustness of the measurement.
The total model uncertainty is the quadratic sum of all contributions. Overall, the uncertainty due to the Dalitz plot amplitude model is small compared to the statistical uncertainty and the experimental systematic uncertainty.

V. INTERPRETATION OF THE RESULTS
The statistical significance of the results is determined by a likelihood-ratio approach. The change in 2 ln L is computed when the CP violation parameters are fixed to zero. The experimental systematic uncertainties and the Dalitz plot amplitude model uncertainties are included by convolution of the likelihood curves. The −2∆lnL curves for sin 2β and cos 2β, and β are shown in Fig. 10. When computing −2∆lnL values for sin 2β and cos 2β, the other observable is fixed to the nominal result. The result for sin 2β agrees within 0.7 standard deviations with the world average of sin 2β = 0.691 ± 0.017 [17] measured fromb →ccs transitions. The measurement excludes the hypothesis of cos 2β ≤ 0 at a p-value of 2.5 × 10 −4 . This corresponds to a significance of 3.7 standard deviations, and thus provides the first evidence for cos 2β > 0. The results exclude the hypothesis of β = 0 • at a p-value of 3.6 × 10 −7 . This corresponds to a significance of 5.1 standard deviations, and thus to an observation of CP violation in B 0 → D ( * ) h 0 decays. The measured value for β is in very good agreement with the preferred solution of the Unitarity Triangle with the world average of (21.9 ± 0.7) • [17]. The second solution of π/2 − β = (68.1 ± 0.7) • is excluded with a pvalue of 2.31 × 10 −13 , corresponding to a significance of 7.3 standard deviations. Therefore, the present measurement reduces an ambiguity in the determination of the parameters of the CKM Unitarity Triangle.

VI. CONCLUSION
In summary, we have measured sin 2β and cos 2β with a time-dependent Dalitz plot analysis of B 0 → D ( * ) h 0 with D → K 0 S π + π − decays. The analysis introduces several improvements over previous related measurements, and new concepts. First, the measurement is performed by a simultaneous analysis of the final data samples collected by the BABAR and Belle experiments, totaling about 1.1 ab −1 and containing about 1 240×10 6 BB pairs collected at the Υ (4S) resonance . The novel combined approach enables the doubling of the statistics available for the measurement, and allows the application of common assumptions and the same D 0 → K 0 S π + π − decay amplitude model simultaneously to the data collected by both experiments. Second, a full Dalitz plot amplitude analysis is performed to derive the D 0 → K 0 S π + π − decay amplitude model directly from a high-statistics e + e − → cc data sample. This enables full control over the model-building process, and the propagation of the D 0 → K 0 S π + π − decay amplitude model uncertainties to those of the CP violation parameters. These approaches lead to improvements in the experimental sensitivity and in the robustness of the measurement.
Moreover, the B 0 → D ( * ) h 0 decays allow a theoretically cleaner determination of the CP -violating phase 2β than the "gold plated" decay modes mediated byb →ccs transitions [84]. Therefore, future more precise measure-ments of B 0 → D ( * ) h 0 decays can provide a new and complementary SM reference for 2β.
The combined BABAR+Belle approach allows the access to an unprecedented large data sample totaling more than 1 ab −1 recorded at c.m. energies of the Υ (4S) resonance and enables a unique experimental precision, in particular, for time-dependent CP violation measurements in the neutral B meson system. Our results underline the importance and discovery potential of future heavy flavor physics experiments operated at high instantaneous luminosity such as the B factory experiment Belle II [32], which is expected to collect a data sample of 1 ab −1 by the year 2020 and is designed to collect 50 ab −1 by the middle of the next decade.